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This  thesis  details  the  system  identification  and  initial  system  validation  of  the  an 
Ultra-Quiet  Vibration  Isolation  Platform  (UQP).  With  the  move  toward  lighter  and  more 
flexible  spacecraft,  the  effects  of  vibration  are  of  immense  concern.  As  natural  or  passive 
damping  becomes  less  effective  in  controlling  undesired  vibrations,  active  vibration 
control  becomes  essential.  The  UQP  uses  a  special  configuration  of  the  six  degree  of 
freedom  Stewart  Platform  with  piezoceramic  strut  actuators  and  geophone  sensors.  This 
combination  gives  an  extremely  sensitive  and  responsive  six  degree-of-freedom  active 
vibration  control  system.  Each  actuator  was  designed  to  be  controlled  independently 
without  coupling  with  other  actuators.  In  order  to  develop  control  laws,  the  plant  must  be 
identified  in  terms  of  system  zeros  and  poles  and  the  uncoupled  design  validated. 
Dynamic  modeling  using  parametric  estimation  methods  can  accurately  describe  a 
complex  system.  Using  parameter  estimation  methods,  models  of  the  actuator  system 
dynamics  were  obtained.  A  simple  lead-lag  controller  was  applied  to  individual  actuators 
then  all  six  actuators  acting  simultaneously  to  verify  system  coupling.  Significant 
interaction  between  base  adjoining  actuators  was  discovered. 
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I.        INTRODUCTION 

A.  PURPOSE  FOR  RESEARCH 

With  increased  use  of  lightweight  materials,  spacecraft  will  become  more  flexible 
and  susceptible  to  excitation  from  sources  of  vibration.  With  loss  of  stiffness  and  rigidity, 
passive  damping  alone  may  be  insufficient  to  attenuate  disturbances  to  an  acceptable 
level.  Active  Vibration  Control  (AVC)  is  rapdily  emerging  as  a  method  of  eliminating  or 
reducing  unwanted  vibration.  The  purpose  of  this  thesis  is  to  provide  a  detailed 
description  and  system  identification  of  the  Ultra  Quiet  Vibration  Isolation  Platform 
(UQP).  It  will  be  the  basis  for  future  research  in  the  area  of  AVC.  The  UQP  is  an 
excellent  sytem  for  AVC  and  can  be  applied  to  a  broad  range  of  disturbances.  A 
disadvantage  of  the  UQP  is  the  complexity  of  the  plant  dynamics  and  kinematics. 
Accordingly  the  complexity  of  the  control  system  design  is  increased.  This  thesis  will 
charaterize  the  plant  for  future  control  design  on  the  UQP. 

B.  SCOPE 

The  "plant"  is  the  UQP  which  consists  of  six  piezoceramic  control  actuators.  This 
thesis  will  use  parameter  estimation  methods  to  create  a  dynamic  model  of  the  UQP 
control  actuators  with  the  ultimate  goal  of  extracting  actuator  zeros,  poles  and  transfer 
functions.  These  actuators  are  designed  to  operate  independently.  Theoretically,  each 
actuator  can  be  controlled  without  coupling  or  interaction  with  the  other  actuators.  This 
assumption  is  of  extreme  importance  as  it  critically  impacts  upon  the  performance  of  the 
system  and  control  system  design.  This  assumption  will  be  checked  for  validity  by 
applying  a  lead-lag  feedback  compensator  to  an  actuator  and  then  to  all  six  actuators 
simultaneously. 


II.       EXPERIMENTAL  SETUP 
A.         BACKGROUND 

In  many  spaceborne  applications  the  dynamics  and  control  of  vibrations  must  be 
addressed  as  a  multiple  degree-of-freedom  (DOF)  problem.  Translations  and  rotations 
about  all  three  axis  must  be  considered.  The  Stewart  Platform  is  the  ideal  mechanism  for 
multiple  DOF  vibration  control  applications. 

1.         The  Stewart  Platform 

The  original  motivation  put  forth  by  Stewart  [Ref.  1]  was  to  design  a 
mechanism  capable  of  simulating  flight  conditions  to  train  pilots.  In  order  to  be  realistic 
it  had  to  be  capable  of  translating  and  rotating  in  three  directions  just  as  a  real  aircraft. 
The  original  configuration  consisted  of  a  triangular  plate  and  a  rigid  parallel  base 
connected  by  six  legs  in  a  "linear  coordinate"  leg  system.  At  each  triangle  vertex  two 
legs  were  attached  in  a  three  DOF  joint  mechanism.  These  legs,  with  hydraulic  actuators, 
were  mounted  to  two  DOF  joints  at  the  base.  The  advantages  of  choosing  this 
configuration  were  inherent  rigidity  and  absence  of  bending  moments.  Additionally  with 
this  configuration  the  only  forces  present  were  in  the  plane  of  the  leg.  A  similiar  six  leg 
arrangement  had  also  been  used  by  a  machine  designed  to  study  tire  to  ground  forces  [Ref 
2].  In  this  system  computer  controlled  jacks  were  used  as  actuators.  A  cubic 
arrangement  was  devised  such  that  the  relationship  of  one  actuator  to  any  other  is  the 
same  (see  Figure  2.1).  This  arrangement  showed  inherent  stability  and  the  capability  for 
accepting  large  moments. 


Figure  2. 1 
Cubic  Arrangement 

Stewart  Platform  mechanisms  have  been  the  subject  of  studies  as  multiple  DOF  parallel 
systems  [Ref.  3  and  Ref.  4:p  46].  The  Stewart  Platform  configuration  has  shown  wide 
applicability  from  motion  simulators  to  robotics  and  now  Active  Vibration  Control.  The 
UQP  employs  the  cubic  configuration  of  the  Stewart  Platform  which  has  an  extremely 
important  property  with  respect  to  AVC.  With  the  exception  of  inertial  loads  and  gravity 
forces,  all  other  forces  are  carried  axially.  The  significance  is  that  if  axial  forces  due  to 
vibration  can  be  eliminated  the  vibration  is  eliminated  [Ref.  4:p.  46]. 


B. 


HARDWARE  CONFIGURATION 


The  experimental  hardware  is  divided  into  two  sections:  the  UQP  and  the  Power, 
Control  and  Analysis  Hardware. 


1.         Ultra-Quiet  Vibration  Isolation  Platform  (UQP) 

The  UQP  is  an  adaptation  of  the  Stewart  Platform  for  AVC  designed  and 
built  by  CSA  Engineering  INC.  It  consists  of  a  circular  plate  attached  to  a  rigid  base  by 
six  variable  length  actuators  (see  Figure  2.2). 


Figure  2.2 
Ultra-Quiet  Vibration  Isolation  Platform 

The  base  consists  of  an  aluminum  plate  supported  by  aluminum  stringers  and  longerons. 
The  base  is  used  to  simulate  a  rigid  spacecraft.  Mounted  inside  one  of  the  "spacecraft" 
bays  is  an  AURA  Bass  Shaker  which  is  used  as  a  known  disturbance  source.  The  entire 
aparatus  is  bolted  to  a  38001b  NEWPORT  optical  table. 


a) 


Control  Actuator  Struts 


The  UQP  uses  six  mutually  orthoganol  struts  to  provide  control 
over  six  degrees  of  freedom.  Each  strut  consists  of  a  piezoceramic  actuator  and  a 
geophone  sensor.   The  use  of  piezoceramics  for  shape  and  vibration  control  is  becoming 


increasingly  widespread.  For  interested  readers,  further  information  on  piezoceramics 
and  their  applications  in  vibration  control  can  be  found  in  Reference  8.  By  taking 
advantage  of  the  cubic  configuration,  all  six  struts  can  be  considered  linear  motion 
actuators  [Ref  1].  Stewart  Platform  based  AVC  mechanisms  using  the  same  cubic 
configuration  with  Terfanol-D  actuators  [Ref.  4]  and  voice  coil  actuators  [Ref.  5]  have 
also  been  developed.  These  platforms  displayed  significant  reduction  of  vibration  over 
passive  means  alone.  Figure  2.3  is  a  basic  diagram  of  each  actuator.  The  actuator  stroke 
is  approximately  50  microns.  In  the  passive  (open  loop)  mode  it  provides  moderate 
damping  to  low  frequency  vibrations.  In  the  active  (closed  loop)  mode  it  provides 
damping  to  higher  frequency  vibrations.  Passive  damping  is  provided  by  a  flexure  with 
damping  material.  An  extremly  sensitive  geophone  measures  velocity.  This  velocity  is 
the  error  signal  fed  into  the  control  law. 


Passive 

Geophone 

Control 

Figure  2.3 
Actuator 


In  general  applications,  geophones  consist  of  wire  coils  supported  by  soft  springs  under 
the  influence  of  a  magnetic  field.  Vibrations  cause  movement  of  the  magnet  relative  to 
the  stationary  coil  inducing  a  voltage  proportional  to  velocity  [Ref.  6:p  449].  The 
stiffness  K  refers  to  the  stiffness  of  the  passive  stage.  The  stiffness  Ka  corresponds  to  the 
piezoceramic  stack  actuator  (PZT). 


b) 


Disturbance  Shaker 


To  create  a  disturbance  source  against  which  the  performance  of 
the  UQP  and  applied  control  can  be  measured,  a  known  disturbance  source  was  mounted 
to  the  base.  A  sinusoidal  waveform  from  a  signal  generator  is  amplified  and  input  to  a 
AURA  Bass  Shaker  (see  Figure  2.4).  This  simulates  a  cyclic  disturbance.  It  has  a 
resonant  frequency  at  42  Hz  which  is  used  as  the  disturbance  frequency  in  initial 
experiments. 


Disturbance 
Shaker 


Figure  2.4 
Disturbance  Shaker 


2.  Power,  Control  and  Analysis 

To  operate  the  UQP  and  extract  meaningful  experimental  data  requires 
several  important  components  (see  Figures  2.5  &  2.6). 


Control  Input-Sensor  Output 


Sensor  Output-  Channel  of  Interest 


PWR  Supply 


AVCS       CK- 


IIHIf-^ 


HP  35665A 


chl _ch2 


KEPCO  Operational 
Amp.  &  Pwr.  Supply 


HP  33 120A  Function/ 
Arb.  Waveform  Generator 


Figure  2.5. 
Power,  Control  and  Analysis  Hardware  Configuration 


Figure  2.6 
High  Voltage  Power  Supply  (HVPS)(top)  and  Active  Vibration  Control  System  (AVCS) 


The  current  configuration  requires  a  MATLAB  capable  PC  with  an  RS232  I/O  port.  The 
control  design  is  done  in  MATLAB.  Once  the  design  is  completed,  the  MATLAB  code  is 
translated  into  C  code.  The  C  code  is  converted  to  machine  language  and  down  loaded 
via  the  I/O  port  to  the  CSA  Active  Vibration  Control  System  (AVCS)  (see  Figure  2.6).  It 
provides  the  interface  for  control  implementation  and  geophone  sensor  output.  Its  main 
component  is  a  digital  signal  processor  (DSP).  Once  the  machine  code  has  been  loaded 
the  AVCS  applies  the  control  to  the  actuators.  It  receives  the  feedback  signal  from  the 
geophones,  processes  and  feeds  it  into  the  control  algorithm.  The  output  is  sent  to  a  HP 
35665A  Dynamic  Signal  Analyzer  for  analysis.  The  CSA  High  Voltage  Power  Supply 
(Figure  2.6)  supports  the  power  equirements  of  the  piezoceramic  actuators  and  is  not 
required  for  use  in  the  open  loop  mode. 
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III.     SYSTEM  IDENTIFICATION 

In  order  to  design  an  effective  compensator  or  controller  for  a  dynamic  system  the 
plant  must  be  well  known  or  characterized.  By  the  use  of  system  identification 
techniques  a  plant  model  can  be  obtained.  The  use  of  a  model  is  extremely  beneficial  in 
cases  when  the  dynamics  of  a  plant  are  not  well  known  or  complex.  Although  articles 
have  been  written  investigating  the  dynamic  and  kinematic  complexities  of  Stewart  type 
platforms  [Ref.  3.],  active  vibration  control  using  Stewart  Platform  based  mechanisms  is 
a  relatively  new  field.  Experimental  dynamic  modeling  can  be  used  to  bypass  the  need 
for  a  complete  theoretic  dynamic  analysis  of  the  UQP.  The  goal  of  identifying  the  UQP 
plant  was  to  obtain  an  accurate  plant  transfer  function  which  can  assist  control  design.  Of 
particular  interest  were  system  modes  or  resonances  and  interaction  or  couplings  among 
actuators.  An  initial  key  assumption  on  the  actuators  is  that  they  each  act  independently 
(uncoupled)  and  can  be  controlled  as  such.  Validation  of  this  assumption  is  one  of  the 
goals  of  identification  of  each  plant 

A.         DYNAMIC  MODELING 

There  are  two  categories  of  system  identification  used  in  dynamic  modeling:  Non- 
parametric  and  Parametric.  Non-parametric  methods  are  used  to  estimate  transient 
response,  spectra  and  frequency  functions  in  order  to  gain  basic  knowledge  of  the  system. 
Parametric  methods  are  used  to  obtain  the  mathematical  parameters  or  elements  used  by 
well  know  system  modeling  algorithms.  These  models  are  broken  down  into  two 
subcatagories.  "Tailor  Made"  models  are  built  based  on  physical  principles  (i.e.  the 
parameters  have  some  type  of  physical  interpretation).  "Ready  Made"  or  "Black  Box" 
models  are  general  in  nature  and  use  parameters  that  have  no  direct  physical 
interpretation.  They  are  used  to  characterize  the  relationship  of  output  to  input  of  a 
system.  By  treating  the  UQP  as  a  "Black  Box"  and  analyzing  the  input-output 
relationship  it  is  possible  to  obtain  an  accurate  model  of  the  system  by  using  a  Ready 
Made  model.  [Ref.  10] 
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1.         Background 

The  concept  behind  the  use  of  transfer  functions  is  to  describe  the 
relationship  between  the  input  and  output  of  a  system  and  can  be  expressed  in  terms  of 
either  continuous  or  discrete  time.  The  use  of  discrete  time  representations  are  preferable 
when  the  data  is  sampled  in  nature. 


u(t) 


Transfer 
Function 


^y(t) 


Figure  3.1 
Input-Output  Relationship 


One  of  the  most  important  tools  used  to  gain  an  understanding  of  the  behavior  of  a  system 
is  by  looking  at  its  impulse  response.  For  a  linear  system  (plant)  the  impulse  response  is 
the  output  when  the  applied  input  is  an  impulse  at  t=0. 


6(t) 


-  y(t)=h(t) 


Figure  3.2 
Impulse  Response  Relationship 

If  the  impulse  response  h(t)  of  a  system  is  known,  the  output  to  a  given  input  can  readily 
be  determined  by  convolution  of  the  input  with  the  impulse  response: 

y(t)  =  u(t)*h(t)  (3.1) 

or: 


y(t)=   I/2(f)««(r-m) 


(3.2) 


The  impulse  response  is  a  key  element  of  non-parametric  system  identification. 

An  essential  element  of  system  identification  is  characterizing  the  output. 
In  an  ideal  system  the  output  would  simply  be  the  input  modified  by  the  system  transfer 
function.  In  physical  systems  however,  the  output  also  contains  elements  resulting  from 
internal  and  external  disturbances  (see  Figure  3.3). 
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Disturbance 


Disturbance 


u(t) 


~d> 


Plant=h(t) 


y(t)=\|/(t)+a(t) 


Where: 


y(t)=Total  Output 
\\f(t)=Undisturbed  Output 
G(i)=Disturbance  term 

Figure  3.3 
In  order  to  develop  a  model  of  the  input  output  relationship,  define: 

y(t)  =  y/(t)  +  cr(t) 

y/(t)  =  G(q,p)u(t) 

a(t)  =  H(q,p)e(t) 
where  £(t)  is  white  noise,  p  is  the  parameter  vector: 


(3.3  a,b,c) 


P  = 


00 


l^J 


(3.4) 


(3.5  a,b) 


and  qr  is  the  shift  operator  based  on  sampling  interval  T: 

qTy(t)  =  y(t+T) 

qjy{t)  =  y(t-T) 
Introduce  the  general  difference  equation: 
y(t)  +  axy{t  -\T)  +  a2y(t  -  27)+. . .+aNy(t  -  NT)  =  b0u(t)  +  b,u{t  -  \T)+. . .+bLu(t  -  LT) 

(3-6) 
where  T  is  the  sampling  interval  of  the  discrete  time  representation.     This  standard 
expression  is  a  specialized  geometric  series  that  relates  input  and  output  of  a  system  using 
past  values  of  input,  output  and  a  set  of  parameters  a;  and  bj.    By  applying  the  shift 
operator  qx,  equation  (3.6)  becomes: 
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XOfl+fltfr1  +a2(lT1  +•  •  -+an<lTn  ]  =  u(t)[b0  +blqTl  +b2qT2  +■  ■  ^n^r"  ] 
y(f)      [^0+^^r1+^2^r2+-"+^^7:"] 

Define: 


(3.7) 


G(^,p)  =  4r  (3-8) 

u(t) 


and  introduce  a  time  delay  of  nk  samples  on  the  input: 


G(q,p)  = 


W7     +b\<lT  +b2<lT  +"-+bn<lT 


(3.9) 


Afo)  [l+ W  +a2qT2  +•  •  •+a/za^rna  ] 

Changing  variables  to  be  consistent  with  the  notation  used  by  Ljung  and  Glad  [ref .  1 1  ] 
and  the  MATLAB  System  Identification  Toolbox  (SITB)  [Ref.  12]: 

u   „-nk  ,l  A—nk—V)  ,  i    A-nk-2),         u      (-nk-nb)"] 
b0qT     +blqKT  +b2qf  J+---+bnqxT 


G(q,P)=B^ 


(3.10) 


Hq)  [l+fiqT1+fl<lT2+-+fnf<lTnf 

Applying  to  equation  (3.3b)  and  equation  (3.6): 

y/(t)  +  fiy/(t-T)+---+fnjr\ff(t-nJT)  =  bou(t-nkT)+---+bnbu(t-(nb  +  nk)T   (3-n) 
Similarly  the  disturbance  term  becomes: 

C(q)-      \l+cml  +c2<lT2  +"-+cnc<lTnC] 
D(q)      [l+dtfT1  +d2qj2  +■  ■  -+dndqjnd  \ 

The  coefficients  bj,  Cj,  d;  and  fi  comprise  the  parameter  vector  p.  These  coefficients 
represent  the  unknown  system  parameters.  Essentially  they  approximate  a  complex 
system's  physical  parameters  but  have  no  direct  correlation  to  any  specific  physical 
quantity.  The  parameters  nb,  nc,  nd  and  nf  characterize  the  order  and  type  of  ready  made 
model.  The  model  constructed  by  equations  (3.3)  through  (3.12)  is  shown  in  equation 
(3.13).  It  represents  the  Box-Jenkins  model1  [Ref.  10]. 

y(t)  =  G(q,p)u(t)  +  H(q,p)£(t)  (3. 1 3) 


1  Named  after  the  statisticians  G.E.P.  Box  and  G.  M.  Jenkins 
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The  following  models  are  variations  of  the  Box-Jenkins  model.  The  most 
simple  one  is  the  case  where  the  disturbance  is  not  modeled. 

y(t)  =  G(q,p)u(t)  +  £(t)  (3.14) 

Equation  (3.14)  is  known  as  the  Output  Error  method.  The  difference  between  the  actual 
output  and  undisturbed  output  is  manifested  in  the  white  noise  term  £(t).  The  next  model 
utilizes  the  same  poles  for  the  input  and  disturbance  (noise)  models  defined  as: 

F(q)  =  D(q)  =  A(q)  =  \  +  aiq-T]+-armq-Tna  (3.15) 

Applied  to  equations  (3.10),  (3.12)  and  (3.13)  yields: 

X0  =  58-tt+gg«W  (316) 

or: 

A(q)y(t)  =  B(q)u(t)  +  C(q)£(t)  (3.17) 

This  is  known  as  the  ARM  AX  (AutoRegression  Moving  Average  eXtra  input)  model. 
The  final  model  that  will  be  described  is  the  AutoRegression  eXtra  input(ARX)  model 
(equation  3.1 8).  [Ref.  10] 

A(q)y(t)  =  B(q)u(t)  +  e(t)  (3.18) 

The  next  step  is  to  describe  how  these  ready  made  techniques  model  a 
system  based  on  past  values  of  output  and  input.  Essentially  these  algorithms  attempt  to 
predict  the  output  based  on  a  given  set  of  input-output  data.  Introduce  the  concept  of 
One-Step- Ahead  prediction  of  output.    From  equations  (3.6)  and  (3. 14),  we  have: 

y(t )  =  -al y(t  -  1)-. .  -aMy(t  -  no)  +  bxu(t  -  nk)+. . ■+bnhu(t  -nk-nb  +  l)  +  e(t)    (3.19) 
Since  £(t)  is  white  noise  and  cannot  be  predicted,  the  prediction  becomes: 

y(t, p)  =  -axy(t - 1)-. .  -a^yit  -  na)  +  bxu{t -  nk)+. . .+^w(r  - nk  -  nb  + 1)       (3.20) 
Expanding  to  the  general  case  of  equation  (3.13)  and  dividing  out  the  noise  transfer 
function: 

H-1  (q,  p)y(t)  =  H~X  (q,  p)G(q,  p)u(t)  +  e{t)  (3.21) 

y(t)  -  y(t)  +  //"'  (q,  p)y(t)  =  H~l  (q,  p)G(q,  p)u(r)  +  e(t)  (3.22) 

y(t )  +  [- 1  +  Hl  (q,  p)]y(t)  =  H~l  (q,  p)G(q,  p)u(t)  +  e(t)  (3.23) 

y(t)  =  [\-H~]  (q,  p)]y(t)  +  H~l  (q,  p)G(q,  p)u(t)  +  e{t)  (3.24) 
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the  prediction  becomes: 

y(t,  p)  =  [1  -  H~l  (q,  p)]y(t)  +  H~x  (q,  p)G(q,  p)u(t)  (3.25) 

Since  equation  (3.25)  only  contains  past  values  of  the  output  y(t)  and  input  u(t)  the 
difference  between  the  prediction  and  the  output  (the  prediction  error)  is: 

E(t,p)  =  y(t)-y(t,p)  =  e(t)  (3.26) 

where  e(t)  is  white  noise.  [Ref.  9:  p.  56  and  Ref.  10:  p.235] 


B.         FREQUENCY  RESPONSE  OF  ACTUATORS 

The  diagram  in  Figure  3.4  is  a  modification  of  the  experimental  set  up  in  Figure 
2.5  used  to  extract  the  input-output  data. 


Control  Input-Sensor  Output 
1 


KEPCO  Operational 
Amp.  &  Pwr.  Supply 


HP  33 120A  Function/ 
Arb.  Waveform  Generator 


PWR  Supply 


AVCS        B- 


ii  ii  in* 


T 


dSPACE 


Excitation  Source 
50  mV  Pink  Noise 

Seijsor  Output-Excited  Channel 


P? 

*?%&&&■•& 


PC 


HP  35665A 
ch^ch^ 


Figure  3.4 
Experimental  Setup 

The  purpose  behind  measuring  the  frequency  response  of  each  actuator  is  to  obtain  both 
basic  insight  into  the  plant  and  a  model  verification  baseline.  Each  actuator  was  excited 
by  a  50  mV  "pink  noise"  source  provided  by  the  HP  Dynamic  Signal  Analyzer  (HPDSA). 
The  "pink  noise"  is  random  noise  which  provides  3dB  roll  off  per  octave.  The  motivation 
to  use  pink  noise  is  to  place  an  equal  amount  of  energy  in  each  octave  band.  The 
response  is  sensed  by  the  geophone  of  the  excited  channel  (actuator)  and  fed  into  the 
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HPDSA.    The  HPDSA  computes  and  plots  the  time  averaged  frequency  response  which 
are  shown  in  Figures  3.5  through  3.10. 
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Figure  3.5 
Frequency  Response  of  Actuator  #1 
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Figure  3.6 
Frequency  Response  of  Actuator  #2 
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Figure  3.7 
Frequency  Response  of  Actuator  #3 
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Figure  3.8 
Frequency  Response  of  Actuator  #4 
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Figure  3.9 
Frequency  Response  of  Actuator  #5 
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Figure  3.10 
Frequency  Response  of  Actuator  #6 


19 


The  plots  show  that  the  plant  transfer  functions  for  the  actuators  are  basically  the  same. 
The  resonant  dynamics  of  the  struts  are  readily  discernible,  specifically  the  actuator 
natural  frequency  at  1 .4  kHz. 


C.         MODEL  STRUCTURE  SELECTION 

The  type  of  system  and  method  of  data  collection  are  major  factors  determining  a 
suitable  model  type  for  a  given  application.  For  example,  the  Output  Error,  OE,  would  be 
best  applied  when  the  disturbances  (which  are  not  modeled)  are  small.  The  ARX  is  the 
most  widely  used  because  it  is  the  simplest  and  serves  as  a  good  baseline  to  select  a 
model.  However,  it  uses  the  same  poles  to  model  noise  as  those  used  to  model  the 
system.  This  can  create  problems  in  system  identification  requiring  the  use  of  higher 
order  models  than  would  otherwise  be  necessary.  The  Box-Jenkins  method  models  the 
system  dynamics  and  disturbances  separately  with  no  common  parameters.  The  ARM  AX 
also  models  the  disturbance  but  uses  the  same  poles  as  the  system  dynamic  model.  It 
differs  from  the  ARX  by  the  addition  of  the  C(q)  polynomial.  The  Box- Jenkins  and 
ARMAX  both  provide  detailed  system  models  and  can  be  extremely  accurate.  The 
decision  on  which  to  use  is  based  on  the  nature  of  the  disturbance.  If  "noise"  enters  the 
process  in  the  early  stages  or  is  carried  through  the  plant  (see  Figure  3.3)  the  ARMAX  is 
preferable  over  the  Box-Jenkins  model  which  is  more  adept  to  characterizing  noises  that 
enter  later  in  the  system.  In  the  case  of  the  UQP,  the  ARMAX  would  appear  to  be  the 
preferred  method  if  there  is  interaction  between  actuators  because  the  resulting 
disturbances  would  enter  into  the  process  early  on.  Additionally  the  disturbance  caused  by 
another  actuator  should  show  similar  dynamics  as  the  plant  under  experimentation.  The 
flow  chart  in  Figure  3.11  represents  the  process  used  in  building  a  model  of  each  actuator. 
As  Figure  3.1 1  shows,  the  process  of  system  identification  and  modeling  is  a  step-by-step 
iterative  process.  The  basic  idea  was  to  develop  and  validate  a  basic  ARX  model  structure 
(i.e.  na,  nb  and  nk)  for  all  six  actuators  and  then  apply  it  to  the  more  complex  ARMAX 
and  Box-Jenkins  methods  and  determine  the  best  model.  The  application  of  this  process 
on  actuator  number  one  will  be  detailed  in  the  following  sections.  [Ref.  101 
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Figure  3.11 
Model  Selection  Flow  Chart 

The  computer  program  used  to  accomplish  this  can  be  found  in  appendix  A. 


1. 


Data  Collection 


The  process  of  system  identification  and  dynamic  modeling  relies  on  a  set 
of  plant  input-output  data.  For  the  UQP  each  actuator  was  provided  an  excitation  and  the 
response  of  all  six  actuators  (as  measured  by  each  respective  geophone)  recorded  (see 
Figure  3.4).  Using  the  HPDSA  as  the  excitation  source,  50mVRMS  "Pink  Noise"  was 
applied  to  the  actuator  via  the  25  pin  connection  on  the  front  of  the  AVCS  (see  Figures 
2.6  and  3.3).  Under  normal  operation  the  ribbon  cable  on  the  front  of  the  AVCS  is 
connected  to  close  the  control  loop.  By  disconnecting  the  cable,  each  actuator  can  be 
accessed  individually.  The  geophone  sensor  output  was  connected  via  a  BNC  interface 
box  to  a  dSPACE  system.  The  dSPACE  was  used  to  gather,  display  and  save  all  six 
channels  of  data  simultaneously.  For  comparison  purposes  the  active  channel  sensor 
output  (actuator  under  excitation)  and  excitation  source  was  also  input  to  the  HPDSA. 
For  each  actuator  the  output  corresponding  to  the  excitation  source  was  measured  at  a 
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sampling  frequency  of  10  kHz  for  two  seconds.  This  gave  20,000  input/output  sample 
points.  The  choice  of  sampling  frequency  was  based  on  the  desire  to  obtain  a  sufficient 
number  of  sample  points  from  which  to  construct  and  validate  a  model.  The  first  10,000 
(input  and  output)  sample  points  are  used  by  the  parameter  estimation  algorithms  to  build 
the  model.  The  last  10,000  points  (of  input)  are  used  in  model  validation.  If  the  sample 
points  used  to  build  the  model  are  also  used  in  validation,  the  true  accuracy  and  validity 
of  the  model  will  be  corrupted. 

a)         Data  Preparation 

Prior  to  using  the  data  to  build  a  model  it  is  necessary  to  pre-treat 
the  data.  The  most  common  factors  adversely  affecting  collected  data  are: 

1 .  Drift,  offset,  trends  and  low  frequency  and/or  periodic  disturbances. 

2.  Outliers  or  faulty  data  points. 

3.  High  frequency  disturbances. 

Before  treating  for  these  possible  problems  the  data  is  given  a  preliminary  analysis  to  see 
if  they  do  in  fact  exist  or  could  be  the  possible  source  of  erroneous  results  [Ref.  9:p  386]. 

External  sources  are  the  main  cause  of  drifts,  trends  and  low 
frequency  disturbances  that  are  either  impossible  or  undesirable  to  model.  This  can  be 
remedied  by  removing  external  disturbances,  using  the  noise  model  to  account  for  the 
disturbances  and/or  using  an  algorithm  to  de-trend  the  data.  Due  to  the  extreme 
sensitivity  of  the  geophone  sensors  it  is  impossible  to  remove  all  the  external 
disturbances.  The  latter  two  options  were  selected.  It  is  assumed  that  the  external 
disturbances  received  by  the  geophones  are  relatively  small  compared  to  those  entering 
into  the  process  early  on  and  will  be  adequately  addressed  by  the  noise  model.  The  offset 
problem  is  a  result  of  the  failure  of  input-output  data  to  correspond  in  a  consistent 
manner.  This  causes  the  model  to  waste  parameters  in  attempt  to  adjust  these  levels. 
This  is  compensated  by  a  MATLAB  function  included  in  the  computer  code  (see 
appendix  A.).  [Ref.  12:pl-63] 
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In  any  data  collection  effort  there  are  obviously  erroneous  values. 
However,  removal  of  these  "outliers"  requires  caution.  Removal  can  either  be  manual 
(i.e.  selecting  bad  values  by  hand)  or  by  some  type  of  automatic  algorithm.  The  best  way 
to  determine  if  outliers  are  problematic  is  to  check  the  model  residuals  for  excessively 
large  values  that  stand  out.  This  method  was  chosen  due  to  the  vast  amounts  of  sampled 
data  used.  [Ref.  12:pl-63] 

The  final  source  of  data  deficiency  is  the  presence  of  high 
frequency  disturbances  outside  the  region  of  interest.  In  this  case  the  primary  concern  is 
the  performance  of  the  actuators  up  to  and  including  the  resonance  frequency.  A  10th 
order  Butterworth  filter  with  a  pass  band  below  approximately  1.5kHz  was  used  to 
eliminate  disturbances  outside  the  range  of  interest  [Ref.  1 1:  p  272]. 

2.  Delay  and  Parameter  Number/Structure  Selection 

Prior  to  applying  one  of  the  aforementioned  parametric  models  to  the 
gathered  data,  the  input  delay  nk,  and  the  structural  elements  na,  nb,  nc,  nd  and  nf  must 
be  determined.  This  was  done  using  the  Model  Structure  Selection  functions  in  the 
MATLAB  SITB.  The  SITB  has  functions  to  calculate  ARX-based  and  OE-based 
structures.  Figure  3.12  represents  a  flow  chart  of  the  algorithm  used  to  select  the  delay 
and  structure.  Since  the  assumption  was  that  the  ARM  AX  model  would  provide  the  best 
model  for  the  UQP,  an  ARX-based  structure  selection  algorithm  was  applied  (see 
appendix  A). 
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Figure  3.12 
Structure  Selection  Algorithm 

The  process  of  selecting  na  and  nb  as  well  as  the  number  of  total  parameters  from  which 
they  are  derived  is  an  iterative  process.  Since  the  input  delay  should  be  the  same  no 
matter  what  the  order  of  the  model  is,  a  second  order  ARX  model  will  serve  as  the 
foundation  of  this  process.  Table  3.1  shows  the  results  of  nk  values  ranging  from  1  to  10 
applied  to  a  second  order  ARX  model. 
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Act. 

nk=l 

nk=2 

nk=3 

nk=4 

nk=5 

nk=6 

nk=7 

nk=8 

nk=9 

nk=10 

1 

-5.5975 

-5.6122 

-5.4704 

-5.4708 

-5.4659 

-5.4686 

-5.4755 

-5.5046 

-5.5184 

-5.5077 

2 

-5.3380 

-5.3175 

-5.2293 

-5.2301 

-5.2297 

-5.2375 

-5.2463 

-5.2602 

-5.2596 

-5.2426 

3 

-4.6214 

-4.4221 

-4.4132 

-4.4391 

-4.4432 

-4.4293 

-4.4357 

-4.4282 

-4.4252 

-4.4256 

4 

-5.5237 

-5.5711 

-5.3791 

-5.3872 

-5.3928 

-5.3946 

-5.4007 

-5.4106 

-5.4054 

-5.3905 

5 

-5.0731 

-5.1164 

-5.0082 

-5.0077 

-5.0044 

-5.0028 

-5.0089 

-5.0239 

-5.0259 

-5.0130 

6 

-5.8464 

-5.8072 

-5.7362 

-5.7301 

-5.7204 

-5.7390 

-5.7443 

-5.7598 

-5.7747 

-5.7610 

Table  3.1 
Input  Delays 

The  numbers  represent  the  logarithms  of  a  quadratic  loss  function  based  on  the  selected 
nk  for  na=nb=2  (second  order  model).  Once  an  input  delay  has  been  selected  a  plot  of 
the  number  of  parameters  used  versus  loss  function  is  produced.  Based  on  the  results  in 
Table  3.1,  the  best  choice  for  an  input  delay  for  actuator  number  one  appears  to  be  nk=2. 
Figure  3.13.  is  the  resulting  plot. 
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Figure  3.13 
Parameters  Vs.  Loss  Function  for  Actuator  #1  nk=2 

The  multiple  results  for  a  given  number  of  parameters  is  due  to  the  fact  that  there  can  be 
several  different  structures  (i.e.  combinations  of  na  and  nb).  The  SITB  has  functions  that 
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automatically  compute  the  structures  based  on  the  minimization  of  Akaike's  Information 
Theoretic  Criterion  (AIC)  and  Rissanen's  Minimum  Description  Length  (MDL)[Ref.  11]. 
After  a  structure  has  been  selected  it  can  be  applied  to  a  parametric  estimation  method  to 
develop  a  model.  Even  with  the  availability  of  functions  capable  of  selecting  the 
optimum  number  of  parameters  based  on  minimization  techniques,  it  can  be  seen  in 
Figure  3.13  that  there  is  a  point  of  diminishing  returns  beyond  which  the  addition  of  more 
parameters  is  of  little  benefit.  Further,  the  addition  of  more  parameters  than  necessary 
will  actually  cause  a  loss  in  model  accuracy.  Referring  to  Figure  3.13  the  structures 
chosen  are  listed  in  table  3.2. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

30 

15  15  2 

AIC  Based 

22 

14  8  2 

MDL  Based 

18 

14  4  2 

Automatic 

Table  3.2 
Structures  For  Actuator  #1  Input  Delay  nk=2 

The  notes  column  provides  information  on  how  the  structure  was  obtained.  AIC  and 
MDL  based  structures  were  automatically  generated  along  with  the  plot  in  Figure  3.13. 
The  term  "Automatic"  means  the  number  of  parameters  in  the  first  column  were  manually 
extracted  from  the  plot  in  Figure  3.13  and  entered  in  to  an  SITB  function  that 
"Automatically"  selects  the  structure.  The  term  "Manual"  refers  to  the  complete  structure 
being  manually  selected. 


3.  Model  Structure  Application  and  Evaluation 

After  determining  the  ARX  model  structure  for  actuator  number  one,  a 
model  was  computed  using  the  algorithm  in  appendix  A.  There  are  numerous  methods 
for  measuring  the  accuracy  and  validity  of  a  model.  The  ones  utilized  to  validate  the 
actuator  models  are  described  as  they  were  applied  to  the  ARX  model  of  actuator  number 
one  using  the  structure  [na=15  nb=15  nk=2]. 
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models  are  described  as  they  were  applied  to  the  ARX  model  of  actuator  number  one 
using  the  structure  [na=15  nb=15  nk=2]. 

a)         Frequency  Response 

One  of  the  quickest  ways  to  determine  a  model's  validity  is  to 
compare  the  Bode  plot  of  the  system  to  that  of  the  model. 
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Figure  3  14 
Frequency  Response  of  ARX[15  15  2]  Model  (Actuator  #1) 


27 


The  model's  frequency  response  is  similar  to  that  of  the  actual  system  (see  Figure  3.14). 
The  areas  of  critical  concern  are  the  model  behavior  at  resonance  and  zero  dB  crossing  for 
magnitude  and  -180  degree  crossing  for  phase.  The  phase  is  reasonably  close  however, 
the  magnitude  is  significantly  different.  The  general  shape  and  resonance  peak  are  in 
agreement  but  there  is  a  distinct  10-20dB  difference  between  the  actuator  and  the  model 

b)         Output  Comparison 

In  evaluating  a  model,  as  stated  earlier,  it  is  important  to  use 
different  input  sample  data  to  insure  an  accurate  assessment  of  the  model.  In  modeling  the 
actuators  a  10,000  sample  model  validation  set  was  held  aside  for  this  puppies.  Applying 
the  validation  input  to  both  the  model  and  the  system  and  comparing  the  outputs,  gives  an 
accurate  measure  of  the  model's  validity  and  how  close  it  truly  simulates  the  actual 
system.  Figure  3.15  compares  the  model  output  to  the  actuator  output  in  response  to  the 
same  input. 
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Figure  3.15 
Actual  Vs.  Model  Output  For  ARX  [15   15  2]  (Actuator  #1) 
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c) 


Auto  and  Cross  Correlation  Functions 


Figure  3.16  provides  plots  of  the  auto  correlation  function  of 
residuals  and  cross-correlation  function  between  the  input  and  residuals  from  the  output. 
The  difference  between  the  model  output  and  actual  output  is  called  the  residual. 
Basically,  it  is  what  is  left  unaccounted  for  by  the  model.  Recalling  equation  (3.26): 

E(t,p)  =  y(t)-y(t,p)  =  £(t) 

where  e(t)  represents  the  residuals.  If  e(t)  is  purely  "white  noise"  it  is  independent  of  the 
input.  If  there  is  correlation  between  e(t)  and  u(t),  there  are  elements  in  the  output  from 
the  input  that  are  not  explained  by  the  model. 
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Figure  3.16 
Auto  and  Cross-Correlation  Functions  for  ARX  [15  15  2]  Model  for  Actuator  #1 


The  horizontal  lines  in  Figure  3.16  indicate  a  99%  confidence  level.  If  the  function 
crosses  these  lines,  there  is  a  correlation  between  £(t)  and  u(t)  at  that  point.  From  Figure 
3.16  there  appears  to  be  correlation  between  e(t)  and  u(t-J).  This  indicates  that  the  input 
delay  might  need  to  be  modified  to  one.  The  auto-correlation  function  is  acceptable.  A 
plot  of  residuals  versus  sample  number  for  100  samples  is  given  in  Figure  3.17. 
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Residual  of  ARX  Model  [15  15  2] 
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Figure  3.17 
Residuals  Vs.  Sample  ARX[15  15  2]  For  Actuator  #1 

Over  this  range  of  samples  the  residuals  are  relatively  low.  For  an  average  output 
magnitude  of  approximately  0.18,  the  average  residual  magnitude  is  on  the  order  of  0.03 
or  roughly  17%  of  the  output  magnitude.  Although  this  percentage  is  a  very  rough  way  of 
comparison  and  cannot  be  considered  by  itself  to  determine  model  validity,  it  does  give  a 
basic  idea  of  how  much  is  left  unexplained  by  the  model.  A  much  better  method  is  to  use 
the  mean  square  fit.  [Ref.10] 


d) 


Zero-Pole  Plot 


The  final  evaluation  tool  is  the  zero-pole  plot  (see  Figure  3.18). 
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Zero  Pole  Plot  ARX  Model  [15  15  2]  Model 
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Figure  3  18 
Zero-Pole  Plot  of  ARX  [15  15  2]  For  Actuator  #1 

The  ellipses  correspond  to  confidence  limits  to  three  standard  deviations  [Ref.  12].  This 
plot  shows  that  there  are  numerous  possible  zero-pole  cancellations  indicating  that  a 
model  of  lesser  degree  is  desirable. 

4.  Model  Application  and  Evaluation  (Data  Prefiltered) 

Based  on  analysis  the  ARX  [15  1 5  2]  is  not  a  very  good  model.  Following 
the  flow  chart  the  same  analysis  will  be  repeated  on  data  pretreated  with  a  10th  order 
Butterworth  filter. 

a)         Frequency  Response 

By  prefiltering  the  data  the  correlation  between  the  model 
frequency  response  and  the  actual  output  frequency  response  is  improved.  However, 
there  is  still  an  approximate  10  dB  difference  between  the  two  responses. 
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Figure  3.19 

Frequency  Response  of  ARX[15  15  2]  Using  Prefiltered  Data  (Actuator#l) 

b)         Output  Comparison 

Figure  3  20  shows  an  improvement  in  the  correlation  between  the 
model  and  actual  output. 
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Plot  of  ARX  MODEL  [15  15  2]  Output  vs.  Actual  output  'Data  Prefiltered* 
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Figure  3.20 
Prefiltered  ARX[15  15  2]  Model  Output  Vs.  Actual  Output 


c) 


Auto  Correlation  and  Cross  Correlation  Functions 


As  a  result  of  prefiltering  there  is  a  degradation  of  the  auto  and 
cross  correlation  functions  (see  Figure  3.21).  However,  the  auto-correlation  of  residuals 
is  relatively  small  For  cross  correlation  functions  where  the  crossing  occurs  in  negative 
lag,  this  is  an  indication  that  the  data  was  possibly  collected  during  feed  back  vice  a 
problem  with  the  model  [Ref.  10:p  284] 
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Figure  3.21 
Auto  and  Cross  Correlation  Functions  of  Prefiltered  ARX[15  15  2]  (Actuator  #1) 
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d) 


Residuals 


Comparing  Figure  3.22  with  Figure  3.20,  the  former  residual  plot 
is  misleading.  This  anomaly  occurs  when  large  order  models  (large  number  of 
parameters)  using  prefiltered  data  are  applied  to  the  SITB  function.  This  function  uses 
statistical  means  to  calculate  the  prediction  error  and  the  residuals.  A  simple  plot  of  the 
actual  output  minus  the  model  output  is  shown  in  Figure  3.23. 
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Figure  3.22 
Residuals  of  Prefiltered  ARX[15  15  2]  (Actuator  #1) 
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Figure  3.23 
Actual  Output  Minus  ARX[15  15  2]  Output  (Actuator  #1) 
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The  average  "residual"  magnitude  for  untreated  and  pretreated  data  are  0.0041  and 
0  0032  respectively,  or  2%  and  1.7%  of  the  average  output  magnitude.  As  previously 
stated,  these  are  very  rough  measures  of  the  accuracy  but  there  is  a  clear  improvement  in 
the  model  by  prefiltering  the  data 


e) 


Zero-Pole  Plots 


The    zero-pole    plot    of    the    prefiltered    model     also     shows 
improvement  over  the  unfiltered  model  (see  Figures  3.18  and  3.24) 


Zero  Pole  Plot  ARX  Model  [15  15  2]  Model  'Data  Prefiltered' 
15 


-0.5 


-1.5         -1  -0.5  0  0.5  1  1.5 


Figure  3.24 
Zero-Pole  Plot  For  Prefiltered  ARX  [15  15  2]  (Actuator  #1) 

However,  there  are  still  numerous  possible  zero-pole  cancellations,  indicating  that  use  of  a 
lower  order  model  is  necessary. 

5.  Model  Structure  Acceptance  for  Actuator  Number  One 

The  goal  of  the  model  structure  application  and  evaluation  process  is  to  use 
the  ARX  based  technique  to  find  an  acceptable  model  structure  that  can  be  applied  to 
more  accurate  ARMAX  and  Box-Jenkins  methods.  Knowing  that  prefiltering  data 
provides  a  better  model,  the  process  is  repeated     Through  the  first  iteration  it  was  also 


35 


determined  that  a  delay  change  would  possibly  be  beneficial.    Returning  to  the  delay 
selection  process,  Figure  3.25  plots  the  loss  function  for  a  given  number  of  parameters. 
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Figure  3.25 
Number  of  Parameters  Vs.  Loss  Function  For  nk=\  (Actuator  #1) 

The  structures  evaluated  are  contained  in  table  3.3. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

30 

15  15  1 

AIC/MDL  Based 

24 

14  10  1 

Automatic 

20 

7  13  1 

Automatic 

18 

7  11  1 

Automatic 

15 

69  1 

Automatic 

13 

67  1 

Automatic 

13 

76  1 

Manual 

Table  3.3 
Structures  For  Actuator  #1  Input  Delay  nk=\ 
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After  a  set  of  candidate  structures  was  obtained  they  were  used  to  generate  ARX  models. 
The  process  of  model  selection  and  evaluation  outlined  in  Figure  3.11  and  described  in 
sections  3  and  4  was  applied  to  each  candidate  structure  based  ARX  model.  The  most 
suitable  candidate  found  was  ARX  [na=l  nb=6  nk=l]  (see  Figures  3.26  though  3.32). 
Table  3.3  lists  the  [na=l  nb-6  nk=Y\(oi  [7  6  1])  structure  as  a  manual  selection.  In  the 
first  iteration  the  [7  13  1]  structure  was  selected.  Upon  repeating  the  selection  process  for 
actuator  number  two,  the  [761]  structure  was  generated  automatically  by  the  algorithm 
(see  appendix  A)  based  on  the  choice  of  13  parameters.  The  performance  of  the  resulting 
model  was  judged  superior  and  the  [7  6  1]  structure  was  applied  to  actuator  number  one 
where  it  also  yielded  the  best  ARX  model  (see  Figure  3.32).  However,  two  design  trade- 
offs were  made.  First,  the  auto  and  cross  correlation  functions  (see  Figure  3.28)  showed 
degraded  performance,  indicating  there  are  still  unaccounted  input  elements  in  the  output. 
However,  it  does  not  affect  the  frequency  response  and  other  criteria  used  to  validate  the 
structure.  The  second  trade-off  was  the  zero-pole  plot  (see  Figure  3.31  ).  The  close 
proximity  of  two  poles  and  two  zeros  indicates  that  a  lower  order  model  would  be 
desirable.  However,  this  proved  to  not  be  the  case.  The  use  of  a  lower  order  model 
showed  declining  performance  over  the  ARX  model  using  the  [761]  structure. 
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Figure  3.26 
Frequency  Response  of  Prefiltered  ARX[7  6  1]  (Actuator  #1) 
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Plot  of  ARX  MODEL  [7  6  1]  Output  vs.  Actual  output  *Data  Prefiltered* 
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Figure  3.27 
Actuator  #1  Output  Vs  ARX[7  6  1]  Model 
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Figure  3.28 
Auto  and  Cross  Correlation  Functions  of  ARX[7  6  1]  Model  (Actuator  #1) 
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Residual  of  ARX  Model  [7  6  1]  *Data  Prefiltered* 
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Figure  3.29 
Residual  Vs.  Sample  for  ARX[7  6  1]  (Actuator  #1) 
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Figure  3.30 
Actual  Output  and  ARX[7  6  11  Model  Output  Difference  (Actuator  #1) 
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Figure  3  31 
Zero-Pole  Plot  for  ARX[7  6  1]  Model  (Actuator  #1) 


41 


Strut  #1  Magnitude 


-50 


o>-100 
Q 


-150 


-200 


10° 
Frequency  (Hz) 

(a)  Magnitude 
Strut  #1  Phase 


"  ■""**-■* 

vl 

\ 

y 

\r 

10° 

Frequency  (HZ) 
(b)  Phase 

Figure  3.32 

ARX  Model  Frequency  Response  Comparison 

Blue=Actuator#l,  Red=[7  6  1],  Magenta=[15  15  1],  Black=[7  13  1] 
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By  increasing  the  delay  the  correlation  functions  improve  but  performance  in  all  other 
respects  is  degraded.  The  99%  confidence  limit  may  also  be  too  constrictive.  Confidence 
limits  of  95%  have  been  used  [ref.  10  p.447]  to  constrain  correlation  functions.  These 
factors  governed  the  decision  to  accept  the  degraded  correlation  functions  as  a  design 
trade  off.  In  comparison  with  models  based  on  various  structures  (Figure  3.32)  this 
model  is  acceptable.  Table  3.4  summarizes  numerical  results. 


Measure 

Result 

Average  Output  Magnitude 

0.1870 

Average  Residual  Magnitude 

0.0064 

Percentage  of  Average  Output 

3.4% 

Average  Output  Difference  Magnitude 

0.0024 

Mean  Square  Fit 

0.1562 

Table  3.4 
Prefiltered  ARX[7  6  1]  Numerical  Results  (Actuator  #1) 

This  process  was  repeated  to  find  the  structures  for  an  ARX  based  model  for  actuator 
numbers  two  through  six.  This  process  involved  the  analysis  of  hundreds  of 
validation/verification  plots  and  numerical  results.  Only  the  validation  plots  of  the 
accepted  structure  for  each  respective  actuator  were  included  in  this  thesis  and  are  shown 
in  appendix  B. 

6.  Model  Structure  Acceptance  for  Actuator  Number  Two 

The  process  outlined  in  section  five  was  repeated  for  actuator  number  two. 
The  plots  resulting  from  the  acceptance  process  are  listed  in  appendix  B  (Figures  A.l 
through  A. 7).  Based  on  a  delay  of  nk=\  from  table  3.1  and  resulting  parameter  number 
plot  (see  appendix  B  Figure  A.l)  the  structures  selected  for  evaluation  are  listed  in  table 
3.5. 
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Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

26 

12  14  1 

AIC  Based 

25 

11  14  1 

MDL  Based 

24 

11  13  1 

Automatic 

22 

11  11  1 

Automatic 

18 

991 

Automatic 

15 

961 

Automatic 

13 

761 

Automatic 

12 

661 

Automatic 

Table  3.5 
Evaluated  Model  Structures  for  Actuator  #2 


After  evaluating  the  graphical  and  numerical  results  of  the  structures  listed  in  table  3.5, 
the  [7  6  1]  structure  was  identified  as  the  best  structure.  It  provides  excellent 
performance  for  a  minimal  number  of  zeros  and  poles.  Figure  3.33  compares  the  [7  6  1] 
based  model  with  the  next  best  models. 


44 


Strut  #2  Magnitude 


Q 


10° 

Frequency  (Hz) 
(a)  Magnitude 
Strut  #2  Phase 


50 

0 

-50 

-100 

-150 

v 

*i 

St 

-200 

/'y 

|v  F 

I, 

10' 


Frequency  (Hz) 
(b)  Phase 

Figure  3.33 

ARX  Model  Frequency  Response  Comparison 

Blue- Actuator,  Red=[7  6  I],  Magenta=[l  1  14  1],  Black=[l  1  111] 
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Measure 

Result 

Average  Output  Magnitude 

0.1906 

Average  Residual  Magnitude 

0.0073 

Percentage  of  Average  Output 

3.8% 

Average  Output  Difference  Magnitude 

0.0025 

Mean  Square  Fit 

0.1736 

Table  3.6 
Numerical  Results  for  Actuator  #2 


7.  Model  Structure  Acceptance  for  Actuator  Number  Three 

The  plots  resulting  from  the  acceptance  process  are  listed  in  appendix  B 
(Figures  B.l  through  B.7).  Based  on  a  delay  of  nk=\  from  table  3.1  and  resulting 
parameter  number  plot  (see  appendix  B  Figure  B.l)  the  structures  selected  for  evaluation 
are  listed  in  table  3.7. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

25 

15  10  1 

Automatic 

24 

15  9  1 

ADC  Based 

24 

14  10  1 

Manual 

17 

89  1 

MDL  Based 

15 

87  1 

Automatic 

13 

85  1 

Automatic 

13 

76  1 

Manual 

Table  3.7 
Evaluated  Model  Structures  for  Actuator  #3 

A  comparison  of  the  selected  structure,  [7  6  1],  and  the  next  best  models  is  shown  in 
Figure  3.34.  Of  particular  interest,  the  best  structures  were  those  resulting  from  non- 
filtered  data.  A  possible  explanation  is  the  absence  of  high  frequency  disturbances  which 
filtering  is  designed  to  remove.  Also  of  note  is  the  fact  that  during  excitation  of  actuator 
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number  three  the  noise  emitted  was  distinctly  different  from  the  other  actuators  when 
excited.  This  could  be  the  result  of  slight  manufacturing  differences. 
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Figure  3.34 

ARX  Model  Frequency  Response  Comparison 

Blue=Actuator#3,  Red=[7  6  1],  Magenta=[15  15  1],  Black=[7  13  1] 
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Measure 

Result 

Average  Output  Magnitude 

0.2298 

Average  Residual  Magnitude 

0.0613 

Percentage  of  Average  Output 

26% 

Average  Output  Difference  Magnitude 

0.0037 

Mean  Square  Fit 

0.2061 

Table  3.8 
Numerical  Results  For  Actuator  #3 


8.  Model  Structure  Acceptance  for  Actuator  Number  Four 

The  plots  resulting  from  the  acceptance  process  are  listed  in  appendix  B 
(Figures  C.l  through  C.7).  Based  on  a  delay  of  nk=\  from  table  3.1  and  resulting 
parameter  number  plot  (see  appendix  B  Figure  C.l)  the  structures  selected  for  evaluation 
are  listed  in  table  3.9. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

26 

11  152 

ADC  Based 

20 

8  122 

MDL  Based 

18 

1082 

Automatic 

12 

482 

Automatic 

26 

11  15  1 

ADC  Based 

24 

11  13  1 

MDL  Based 

17 

89  1 

Automatic 

13 

76  1 

Manual 

Table  3.9 

Evaluated  Model  Structures 

Similar  to     the  case  of  actuator  number  one,  the  algorithm  selected  a  delay  of  two. 
However,  the  best  results  were  achieved  by  using  a  lower  delay.    The  penalty  paid  was 
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degraded  auto  and  cross  correlation  functions.  This  was  again  taken  as  a  design  trade  off. 
The  ARX  structure  [na=l  nb=6  nk=\]  again  proved  to  be  the  best  choice  (see  Figure 
3.35).  Numerical  results  are  contained  in  table  3.10. 
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Figure  3.35 

ARX  Model  Frequency  Response  Comparison 

Blue=Actuator#4,  Red=[7  6  1],  Magenta=[ll  13  1],  Black=[ll  111] 
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Measure 

Result 

Average  Output  Magnitude 

0.3052 

Average  Residual  Magnitude 

0.0066 

Percentage  of  Average  Output 

2.1% 

Average  Output  Difference  Magnitude 

0.0036 

Mean  Square  Fit 

0.2662 

Table  3.10 
Numerical  Results  for  Actuator  #4 


9.  Model  Structure  Acceptance  for  Actuator  Number  Five 

The  plots  resulting  from  the  acceptance  process  are  listed  in  appendix  B 
(Figures  D.l  through  D.7).  Based  on  a  delay  of  nk=\  from  table  3.1  and  resulting 
parameter  number  plot  (see  appendix  B  Figure  D.  1 )  the  structures  selected  for  evaluation 
are  listed  in  table  3.12. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

28 

13  15  2 

ADC  Based 

27 

12  15  2 

MDL  Based 

22 

10  122 

Automatic 

18 

1262 

Automatic 

28 

13  15  1 

ADC/MDL  Based 

24 

15  9  1 

Automatic 

15 

87  1 

Automatic 

13 

76  1 

Manual 

Table  3.11 
Evaluated  Model  Structures 


The  validation  of  actuator  number  five  is  virtually  the  same  as  actuator  number  four  (see 
appendix  B  Figures  C.2  through  C.7,  D.2  through  D.7  and  tables  3.10  and  3.12). 
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Figure  3.36 

ARX  Model  Frequency  Response  Comparison 

Blue= Actuator  #5,  Red=[7  6  1],  Magenta=[15  9  1],  Black=[13  15  1] 
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Measure 

Result 

Average  Output  Magnitude 

0.3070 

Average  Residual  Magnitude 

0.0078 

Percentage  of  Average  Output 

2.5% 

Average  Output  Difference  Magnitude 

0.0032 

Mean  Square  Fit 

0.2830 

Table  3.12 
Numerical  Results  for  Actuator  #5 


10.        Model  Structure  Acceptance  for  Actuator  Number  Six 

The  plots  resulting  from  the  acceptance  process  are  listed  in  appendix  B 
(Figures  E.l  through  E.7).  Based  on  a  delay  of  nk=\  from  table  3.1  and  resulting 
parameter  number  plot  (see  appendix  B  Figure  E.  1 )  the  structures  selected  for  evaluation 
are  listed  in  table  3.13. 


Number  of  Parameters 

Structure  [na  nb  nk] 

Notes 

30 

15  15  1 

ADC/MDL  Based 

24 

12  12  1 

Automatic 

20 

11  9  1 

Automatic 

18 

99  1 

Automatic 

13 

76  1 

Manual 

Table  3.13 
Evaluated  Model  Structures 

The  comparison  of  the  best  model  structures  is  shown  in  Figure  3.37  .  The  [7  6  1] 
structure  selected  for  the  other  five  actuators  is  the  best  choice  for  actuator  number  six. 
Numerical  results  are  listed  in  table  3.14. 
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Figure  3.37 

ARX  Model  Frequency  Response  Comparison 

Blue= Actuator  #6,  Red=[7  6  1],  Magenta=[15  15  1],  Black=[12  12  1] 
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Measure 

Result 

Average  Output  Magnitude 

0.1551 

Average  Residual  Magnitude 

0.0055 

Percentage  of  Average  Output 

3.5% 

Average  Output  Difference  Magnitude 

0.0023 

Mean  Square  Fit 

0.1316 

Table  3.14 
Numerical  Results  for  Actuator  #6 
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D.         APPLICATION  OF  SELECTED  STRUCTURE 

The  selection  of  the  same  basic  model  structure  na=l  nb=6  nk=\  for  all  six 
actuators  is  an  indication  that  the  plant  dynamics  of  the  actuators  can  be  assumed  to  be 
identical.  With  the  exception  of  actuator  number  three,  all  received  the  best  results  from 
prefiltered  data.  The  design  trade-offs  were  the  degraded  auto  and  cross  correlation 
functions  and  the  relative  proximity  of  two  poles  and  two  zeros.  This  basic  structure  was 
used  to  develop  more  accurate  ARMAX  and  Box-Jenkins  models.  Although  the  Box- 
Jenkins  is  the  most  accurate  parametric  estimation  method,  it  was  anticipated  that  the 
ARMAX  would  give  the  best  results  based  on  the  nature  of  disturbances  and  on  when 
they  enter  into  the  system  (see  Figure  3.3). 

The  application  and  evaluation  process  is  the  same  as  that  used  to  evaluate 
selected  structures  on  the  ARX  model.  A  second  order  noise  model  was  used  for  both  the 
ARMAX  and  Box-Jenkins  models.  The  standard  choices  are  2nd  and  4th  order  models 
[Ref.  11].  Both  were  applied  and  the  2nd  order  yielded  a  better  result.  After  the 
application  of  the  ARMAX  and  Box-Jenkins  based  models  the  overall  performance  of  the 
ARX  and  ARMAX  models  was  decidedly  better.  The  only  exception  was  on  actuator 
number  two.  This  will  be  addressed  in  section  two.  Based  on  the  results  obtained  thus 
far  indicating  the  actuator  plants  are  essentially  the  same  and  the  excellent  performance  of 
the  ARX  and  ARMAX  models,  only  the  evaluation  plots  for  the  ARMAX  model  are 
included  and  can  be  found  in  appendix  C.  Summary  and  comparison  results  are  listed  for 
each  actuator. 


1.  Actuator  Number  One 

A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.38.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  A.l  through  A.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.15. 
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Figure  3  38 

Model  Frequency  Response  Comparison 

Blue= Actuator  #1,  Red=ARMAX[7  6  2  1],  Magenta=ARX[7  6  1],  Black=BJ[6  2  2  7  1] 
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Measure 

ARX[7  6  1] 
(Data  Prefiltered) 

ARMAX[7  6  2  1] 
(Data  Prefiltered) 

BJ[6227  1J 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.1870 

0.1870 

0.1870 

Ave.  Residual  Magnitude 

0.0064 

0.0036 

0.0094 

Percent  of  Output  Mag. 

3.4% 

2% 

5% 

Ave.  Output  Difference 

0.0024 

0.002 

0.0025 

Mean  Square  Fit 

0.1562 

0.1589 

0.1568 

Table  3.15 
Numerical  Results  for  Actuator  #1 

Based  upon  the  above  results  the  ARX[7  6  1]  model  was  determined  to  best  fit  the 
performance  of  actuator  number  one.  The  choice  between  the  ARX  and  ARMAX  is  a 
close  one.  This  indicates  the  assumption  that  the  disturbance  or  noise  entered  the  system 
early  and  is  accurately  modeled  by  the  dynamics  (poles)  of  the  plant.  The  addition  of  the 
C(q,9)  polynomial  is  not  necessary.  The  improved  performance  of  the  ARX  and 
ARMAX  models  over  the  Box-Jenkins  further  supports  the  assumption  regarding  the 
nature  of  the  disturbance. 


2. 


Actuator  Number  Two 


A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.39.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  B.l  through  B.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.16. 
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Figure  3.39 
Model  Frequency  Response  Comparison 
BIue= Actuator  #2,  Red=ARMAX[7  6  2  1],  Magenta- ARX[7  6  I],  BIack=BJ[6  2  2  7  1] 
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Measure 

ARX[7  6  1] 
(Data  Prefiltered) 

ARMAX[7  6  2  1] 
(Data  Prefiltered) 

BJ[62  27  1] 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.1906 

0.1906 

0.1906 

Ave.  Residual  Magnitude 

0.0073 

0.003 

0.0114 

Percent  of  Output  Mag. 

3.8% 

1.6% 

6% 

Ave.  Output  Difference 

0.0025 

0.0027 

0.0023 

Mean  Square  Fit 

0.1736 

0.1705 

0.1616 

Table  3.16 
Numerical  Results  for  Actuator  #2 


The  Box-Jenkins  model  showed  excellent  frequency  response  correlation,  auto  and  cross 
correlation  functions  and  some  improved  numerical  results.  However,  due  to  the  close 
performance  of  the  ARMAX  and  Box-Jenkins  models,  the  large  uncertainty  ellipses  in 
the  zero-pole  plot  (see  appendix  C)  the  ARMAX  model  was  selected. 


3. 


Actuator  Number  Three 


A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.40.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  C.l  through  C.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.17. 
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Figure  3.40 
Model  Frequency  Response  Comparison 
Blue= Actuator  #3,  Red=ARMAX[7  6  2  1],  Magenta=ARX[7  6  1],  Black=BJ[6  2  2  7  1] 
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Measure 

ARX[7  6  1] 
(Data  not  filtered) 

ARMAX[7  6  2  1] 
(Data  Prefiltered) 

BJ[62  27  1] 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.2298 

0.2298 

0.2298 

Ave.  Residual  Mag. 

0.0613 

0.0040 

0.0215 

Percent  of  Output  Mag. 

26% 

1.7% 

9.3% 

Ave.  Output  Difference 

0.0037 

0.0040 

0.0036 

Mean  Square  Fit 

0.2061 

0.2031 

0.2000 

Table  3.17 
Numerical  Results  for  Actuator  #3 


Although  the  non-filtered  ARX  model  was  chosen  over  the  prefiltered  ARX  model,  it 
appears  that  the  performance  of  the  prefiltered  ARMAX  is  much  better.  This  can  be 
explained  by  possible  manufacturing  variation  in  actuator  number  three  which  requires 
the  use  of  the  C(q,0)  polynomial  to  model  the  disturbance. 


4. 


Actuator  Number  Four 


A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.41.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  D.l  through  D.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.18. 
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Figure  3.41 
Model  Frequency  Response  Comparison 
Blue- Actuator  #4,  Red=ARMAX[7  6  2  1],  Magenta=ARX[7  6  1],  Black=BJ[6  2  2  7  1] 
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Measure 

ARX[7  6  1] 
(Data  Prefiltered) 

ARMAX[7  6  2  1] 
(Data  Not  Prefiltered) 

BJ[6227  1] 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.3052 

0.3052 

0.3052 

Ave.  Residual  Mag. 

0.0066 

0.0335 

0.0099 

Percent  of  Output  Mag. 

2.1% 

11% 

3.2% 

Ave.  Output  Difference 

0.0036 

0.0034 

0.0034 

Mean  Square  Fit 

0.2662 

0.2661 

0.2572 

Table  3.18 
Numerical  Results  for  Actuator  #4 

This  case  appears  similar  to  actuator  number  one  where  the  disturbance  is  accurately 

characterized  by  the  system  pole  model  without  the  need  for  the  C(q,9)  polynomial. 

5.  Actuator  Number  Five 

A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.42.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  E.l  through  E.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.19. 
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Strut  #5  Magnitude 
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Figure  3.42 
Model  Frequency  Response  Comparison 
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Measure 

ARX[7  6  1] 
(Data  filtered) 

ARMAX[7  6  2  1] 
(Data  Prefiltered) 

BJ[6  2  27  1] 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.3070 

0.3070 

0.3070 

Ave.  Residual  Mag. 

0.0078 

0.0035 

0.0151 

Percent  of  Output  Mag. 

2.5% 

1.1% 

4.9% 

Ave.  Output  Difference 

0.0032 

0.0033 

0.0034 

Mean  Square  Fit 

0.2830 

0.2688 

0.2698 

Table  3.19 
Numerical  Results  for  Actuator  #5 


6. 


Actuator  Number  Six 


A  graphical  comparison  of  ARX,  ARMAX  and  Box-Jenkins  models  is 
shown  in  Figure  3.43.  Validation  plots  for  the  ARMAX  [7  6  2  1]  model  are  included  as 
Figures  F.l  through  F.6  in  appendix  C.  Numerical  results  are  contained  in  table  3.20. 
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Figure  3.43 
Model  Frequency  Response  Comparison 
Blue= Actuator  #6,  Red=ARMAX[7  6  2  1],  Magenta=ARX[7  6  1],  Black=BJ[6  2  2  7  1] 
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Measure 

ARX[7  6  1] 
(Data  not  filtered) 

ARMAX[7  6  2  1] 
(Data  Prefiltered) 

BJ[62  2  7  1] 
(Data  Prefiltered) 

Ave.  Output  Magnitude 

0.1551 

0.1551 

0.1551 

Ave.  Residual  Mag. 

0.0055 

0.0035 

0.0086 

Percent  of  Output  Mag. 

3.5% 

2.2% 

5.5% 

Ave.  Output  Difference 

0.0023 

0.0020 

0.0019 

Mean  Square  Fit 

0.1316 

0.1312 

0.1356 

Table  3.20 
Numerical  Results  for  Actuator  #6 


E.         SUMMARY  OF  RESULTS 


Table  3.21  presents  a  summary  of  the  model  selection  process. 


Actuator 

1 

2 

3 

4 

5 

6 

Model 
[Structure] 

ARX 

[7  6  1] 

ARM  AX 

[7  6  2  1] 

ARM  AX 

[7  6  2  1] 

ARX 

[7  6  1] 

ARM  AX 
[7  6  2  1] 

ARM  AX 
[7  6  2  1] 

Ave.  Out.  Mag. 

0.1870 

0.1906 

0.2298 

0.3052 

0.3070 

0.1551 

Ave.  Res.  Mag. 

0.0064 

0.003 

0.0040 

0.0066 

0.0035 

0.0035 

%  of  Out. 
Mag. 

3.4% 

1.6% 

1.7% 

2.1% 

1.1% 

2.2% 

Ave.  Out. 
Dive. 

0.0024 

0.0027 

0.0040 

0.0036 

0.0033 

0.0020 

Mean  Sq.  Fit 

0.1562 

0.1705 

0.2031 

0.2662 

0.2688 

0.1312 

Table  3.21 
Summary  of  Numerical  Results 

F.         POLES  AND  ZEROS 

The  ultimate  goal  in  identifying  the  plant  of  the  platform  is  to  obtain  the  poles  and 
zeros  of  the  plant  transfer  function  which  can  be  used  in  design  of  an  active  vibration 
control  law  for  the  platform.  Tables  3.22  and  3.23  present  the  poles  and  zeros  extracted 
from  the  parameter  estimation  algorithm  listed  in  appendix  A. 
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1 

2 

3 

0.1252+0.9199i 

0.1080+0.9667i 

0.1198+0.9204i 

0. 1252-0.9 199i 

0.1080-0.9667i 

0.1198-0.9204i 

1.0137 

1.0222 

1.0204 

0.5371+0.6167i 

0.5887+0.6486i 

0.4925+0.6206i 

0.5371  -0.6 167i 

0.5887-0.6486i 

0.4925-0.6206i 

(a)  Actuators  #1  through  #3 

4 

5 

6 

0.0674  +  0.6235i 

0.1971+0.8593i 

0.228  l+0.7945i 

0.0674  +  0.6235i 

0.1971-0.8593i 

0.228  l-0.7945i 

-1.2095 

1.0098 

1.0032 

1.0010 

-0.5566 

-0.1319+0.6080i 

-0.5568 

0.2913 

-0.1319-0.6080i 

(b)  Actuators  #4  through  #6 

Table  3.22 
Actuator  Transfer  Function  Poles 
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1 

2 

3 

0.1661  +0.9362i 

0.1122  +  0.9803i 

0.1410 +  0.9402i 

0.1661  -0.9362i 

0.1122-0.9803i 

0.1410 +  0.9402i 

0.8743 

0.9035 

0.7231 

0.5932 +  0.759H 

0.5645  +  0.810H 

0.5584+0.7771i 

0.5932  -0.759H 

0.5645- 0.8 lOli 

0.5584  -0.777 li 

0.7101  +0.5677i 

0.7220  +  0.6035i 

0.5562  +  0.5203i 

0.7101  +0.5677i 

0.7220  -  0.6035i 

0.5562  +  0.5203i 

(a)  Actuators  #1  through  #3 

4 

5 

6 

0.6905  +  0.6666i 

0.1657 +  0.9000i 

0.1138  +  0.8029i 

0.6905  -  0.6666i 

0.1657 -0.9000i 

0.1138-0.8029i 

0.6386 

0.8267 

0.6405 

-0.3437  +  0.4468i 

0.6710 +  0.6870i 

0.5522  +  0.7348i 

-0.3437  -  0.4468i 

0.6710  -0.6870i 

0.5522  -0.7348i 

0.2640  +  0.453  li 

0.4532  +  0.4747i 

0.3908  +  0.4466i 

0.2640- 0.453  li 

0.4532  -  0.4747i 

0.3908  -0.4466i 

(b)  Actuators  #4  through  #6 

Table  3.23 
Actuator  Transfer  Function  Zeros 


A  review  of  the  poles  and  zeros  shows  that  zeros  are  very  similar,  especially  among 
actuators  one,  two,  three  and  six.  The  magnitude  of  the  real  parts  of  zeros  for  actuators 
four  and  five  appear  to  vary  more  than  the  others.  This  is  also  true  for  the  poles.  This 
was  unexpected  and  no  explanation  is  offered.  Any  expected  difference  would  have  been 
attributed  to  actuator  number  three  based  on  previous  results  and  analysis.  However, 
actuator  number  three  is  in  relative  concurrence  with  the  other  actuators. 
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G.        TRANSFER  FUNCTIONS 

The  actuator  transfer  functions  are  as  follows; 

ACTUATOR  Number  One 

13.536256  -3 1.64935s  +42.53  1754  -39.488653  +  22.801252  -7.90855 
57  -3.813256  +  7.77935s  -10.353254  +  9.7 10253  -6.390752  +  2.73365-0.6064 

ACTUATOR  Number  Two 

1 5.1 70956  -  36.64645s  +  5 1.459654  -  49.929653  +  30.859952  - 1 1.25845 
57  -3.701056  +  7.56925s  -10.251 154  +  9.948253  -6.876252  +3.13325-0.7594 

ACTUATOR  Number  Three 

2 1.687556  -  48.690455  +  64.5 1 9854  -  59.848553  +  33.837052  - 1 1.96735 
57  -  3.234256  +  6.067955  -  7.502254  +  6.638053  -  4.1 17552  + 1.64545  -  0.3361 

ACTUATOR  Number  Four 

13.487556  -30.799355  +40.403254  -36.90 1753  +2 1.223652  -7.54895 
57  -3.927056  +8.095355  -10.807654  +  10.133853  -6.634852  +2.80155-0.6022 

ACTUATOR  Number  Five 

8.5 1 3 156  +  3.982355  -  9.423754  -  4.224253  +  0.9300252  -  0.1 9505 
57  -  2.035956  + 1.74355  -  0.482454  +  0.0 1 7553  -  0.277352  +  0.25775  -  0.1 089 

ACTUATOR  Number  Six 

9.676856  - 1 1.56925s  + 1 1.059 154  -  9.1 85853  +  2.522752  -  2.567 15 
57  -2Z754256  +4.501 155  -4.774454  +3.658253  - 1.969252  +0.681 15-0.1253 

There  is  greater  concurrence  among  actuators  one,  two,  three  and  six.  However,  overall 
comparison  of  zeros,  poles  and  transfer  functions  indicates  that  all  six  plants  are  very 
similar. 
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IV.      ACTUATOR  COUPLING  AND  VIBRATION  CONTROL 

APPLICATION 


A. 


ACTUATOR  COUPLING 


With  the  dynamic  models  constructed  in  Chapter  Three  it  is  now  possible  to  model 
and  investigate  the  interaction  between  actuators  to  see  if  indeed  they  can  be  controlled 
independently  The  basic  ARX  method  has  been  proved  to  be  quite  accurate  in  modeling 
the  plant  transfer  function  of  each  actuator  To  determine  the  presence  of  actuator 
coupling  each  actuator  was  excited  with  50mV  pink  noise  and  the  output  response  of  all 
six  actuators  measured.  The  input-output  data  gathered  was  introduced  to  the  ARX[7  6 
1]  model  and  the  resulting  model  frequency  response  was  plotted.  Figures  4  1  through  4.5 
show  the  results  from  the  excitation  of  actuator  number  one. 
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Figure  4  1 
Frequency  Response  Actuator  #1  Vs.  Actuator  #2 
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Figure  4.2 
Frequency  Response  Actuator  #  1  Vs.  Actuator  #3 
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Figure  4.3 
Frequency  Response  Actuator  #1  Vs  Actuator  #4 
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Figure  4.4 
Frequency  Response  Actuator  #1  Vs.  Actuator  #5 
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Figure  4.5 
Frequency  Response  Actuator  #1  Vs.  Actuator  #6 
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In  Figures  4. 1  through  4.4  there  is  discernible  correlation  between  the  frequency 
response  of  actuator  number  one  and  the  frequency  response  of  actuators  two  through  five 
to  the  excitation  of  actuator  number  one.  Of  particular  interest  is  the  response  at  the 
resonance  frequency  of  1.4  kHz.  On  average  there  is  approximately  a  10  dB  difference  at 
1.4  kHz  for  these  actuators.  In  the  frequencies  below  1kHz  the  response  is  in  the 
negative  dB  range.  This  would  indicate  that  the  passive  stage  is  performing  well  in 
eliminating  the  low  frequency  disturbances.  Although  coupling  is  present,  it  is  relatively 
low  in  comparison  to  the  response  of  actuator  one.  However  this  is  not  the  case  with 
actuator  number  six.  As  in  actuators  two  through  five  the  response  below  1  k  is  small. 
Above  1kHz,  particularly  at  1.4  kHz  the  response  is  close  to  a  mirror  image  of  the 
response  of  actuator  number  one.  This  suggests  that  there  is  significant  high  frequency 
coupling 

B.         CONTROL  APPLICATION  AND  COUPLING  VERIFICATION 

Given  the  results  of  the  previous  section  it  is  necessary  to  verify  on  the  UQP. 
This  was  accomplished  by  applying  a  control  law  to  the  platform  and  checking  for  any 
evidence  of  coupling.  The  first  evaluation  conducted  consisted  of  controlling  an 
individual  actuator  and  the  second  applied  control  to  all  six  actuators  simultaneously. 
There  are  numerous  control  candidates  for  use  in  AVC.  For  random  broad  band 
noise/vibrations,  feedback  control  is  preferable.  Feed  forward  controllers  are  good 
candidates  when  a  specific  frequency  of  the  disturbance  is  known.  Local  (decoupled) 
force  feedback  [Ref.  5],  adaptive  vibration  control  using  two  least  mean  square  filters 
[Ref.  4],  and  absolute  velocity  feedback  [Ref.  6]  are  examples  that  have  been  used  in 
AVC.  The  latter  of  the  three  also  used  a  geophone  for  sensing.  The  control  law  currently 
programmed  for  the  system  is  a  lead-lag  compensator.  This  type  of  compensator  is  a 
variant  of  derivative  and  integral  control.  Lead  compensation  has  the  effect  of  lowering 
rise  time  and  decreasing  overshoot  while  lag  compensation  is  used  mainly  to  gain  steady 
state  accuracy  of  the  system[Ref.  7].  A  lead-lag  controller  was  provided  by  CSA 
Engineering  and  was  used  to  ensure  the  operation  of  each  individual  strut.   This  control 
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law  was  modified  to  apply  the  same  lead-lag  compensator  to  all  six  actuators 
simultaneously  based  on  the  initial  assumption  of  independently  controllable  actuators. 
The  results  previously  obtained  indicate  that  this  cannot  be  done. 

1.  Compensator  Transfer  Function 

A  Bode  Plot  of  the  magnitude  and  phase  of  the  compensator  (See  Figure 
4.6)  shows  an  infinite  gain  margin  and  a  phase  margin  of  approximately  60  degrees. 
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Compensator  Bode  Plot 
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A  bode  plot  of  the  open  loop  tranfer  function  is  provided  in  Figure  4.7. 
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Figure  4.7 
Open  Loop  Transfer  Function  (Compensator* Actuator) 

Finally,  a  magnitude  plot  in  terms  of  dB  of  closed  to  open  loop  is  shown  in  Figure  4.8. 
From  this  plot  it  is  observed  that  the  best  performance  for  this  compensator  will  be 
between  12  Hz  to  48  Hz  with  maximum  reduction  of  32dB  at  25  Hz.  A  reduction  of 
approximately  25  dB  can  be  seen  at  42  Hz.  This  controller  was  designed  with  an 
approximate  bandwidth  of  100  Hz.  The  spike  at  approximately  10  Hz  is  a  faulty  value 
resulting  from  system  noise  or  an  error  in  data  collection. 
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Figure  4.8 
Closed  Loop/Open  Loop 

2.  Control  Application 

The  controller  described  in  section  one  was  first  applied  to  actuator 
number  two  and  then  the  .modified  controller  was  applied  to  all  six  actuators 
simultaneously.  Table  4.1  details  the  equipment  preparation  and  settings. 
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Equipment 

Setting 

HP  33120A  Function  Generator 

Waveform                            Sinusoid 
Frequency                            42  Hz 
Amplitude                            500mv  peak- 
to-peak 

KEPCO  Bipolar  Operational  Power 
Supply/Amplifier 

Current  Control  Mode               +/- 1  Amp 

Table  4.1 
Equipment  Settings 


3.         Experiment 


The  first  task  is  to  look  at  the  platform  in  open  loop  mode.  From  this  the 
local  environment  can  be  assessed.  The  plots  in  Figures  4.9  and  4.10  show  what  is 
sensed  by  the  geophone  on  actuator  two  with  the  42  Hz  disturbance  and  no  control  (open 
loop). 
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Figure  4.9 
Ambient  PSD  (0-5kHz) 
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Open  Loop-42Hz  Disturbance 
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Figure  4. 10 
Ambient  PSD  (0-200Hz) 

Of  note  in  Figure  4.9  is  the  1.4  kHz  resonance  frequency.  This  is  the  result  of  power 
being  applied  to  the  actuator  from  the  HVPS.  At  approximately  30  Hz  there  is  an 
environmental  disturbance  of  unknown  origin.  In  the  open  loop  mode  the  42Hz 
disturbance  is  clearly  visible  (see  Figure  4.10).  Figure  4.11  shows  the  power  spectral 
density  (PSD)  of  the  sensor  output  for  actuator  number  two  for  open  and  closed  loop 
operations  (both  single  and  all  six  actuators).  With  only  power  applied  to  the 
piezoceramic  actuator  (no  control),  excitation  of  the  resonance  frequency  is  visible  as  the 
blue  trace  in  Figure  4.11(a).  Closing  the  control  loop  for  single  and  multi-actuator 
control  increases  the  excitation.  Figure  4.12  shows  the  sensor  output  for  actuators  two, 
three,  and  six  with  control  applied  to  actuator  number  two.  Actuator  number  three  shares 
the  same  base  node  with  actuator  number  two.  The  coupling  between  actuators  two  and 
three  can  be  seen  in  Figure  4.12(a).  The  application  of  control  excites  the  resonance  of 
all  three  actuators,  but  the  interaction  between  actuators  two  and  three  is  the  most 
pronounced. 
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(a)  PSD  (0-5kHz)  Sensed  By  Geophone  on  Actuator  #2 

Blue-Open  Loop;  Magenta-Strut  #2  Closed  Loop;  Red-All  6  Struts  Closed  Loop 
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(b)  PSD  (0-200Hz)  Sensed  By  Geophone  on  Actuator  #2 

Figure  4.11 
Open  Vs.  Closed  Loop  Performance  PSD 


82 


Blue-#2  Geophone;  Magenta-#3  Geophone;  Red-#  6  Geophone 
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(a)  PSD  (0-5kHz)  Sensor  Outputs 
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(b)  PSD  (0-200Hz)  Sensor  Outputs 

Figure  4.12 
Actuator  #2  Closed  Loop  Performance  PSD 
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V.       SUMMARY 

A.         CONCLUSIONS 

The  process  of  system  identification  provided  valuable  insight  into  the  UQP.  By 
treating  the  UQP  as  a  "Black  Box"  and  using  ready-made  models  the  complex  system 
dynamics  were  overcome  and  very  accurate  models  of  the  system  were  developed.  Some 
possible  reasons  for  differences  between  the  models  and  the  actual  system  can  be 
explained  by  non-removal  of  faulty  values.  Due  to  the  number  of  samples  used  to  build 
and  validate  the  models,  hand  smoothing  of  the  data  was  not  practical.  Based  on  the 
average  residuals  for  each  actuator  being  less  than  4%  of  the  average  model  output, 
outliers  did  not  pose  a  significant  problem.  In  comparing  the  validation  data  (numerical 
and  graphical),  actual  frequency  responses,  and  zeros,  poles  and  transfer  functions  there 
are  distinct  similarities.  The  fact  that  all  six  are  described  by  seven  poles  and  six  zeros 
was  anticipated.  However,  some  variation  among  actuators  was  expected.  There  is 
strong  correlation  of  actuators  one,  two,  three  and  six  are.  Although  the  numerical  values 
for  the  zeros,  poles  and  transfer  function  of  actuators  four  and  five  did  not  correlate  with 
the  other  actuators  as  strongly,  they  still  retained  the  same  pole-zero  structure.  Based  on 
this  it  can  be  stated  that  for  the  purposes  of  characterizing  the  system  and  control  design, 
all  six  actuators  are  identical. 

The  sensitivity  of  the  geophones  resulting  in  the  introduction  of  external 
disturbances  could  be  a  source  of  inaccuracy,  however,  they  are  considered  small  in 
comparison  to  the  disturbances  generated  by  the  other  actuators.  The  emergence  of  ARX 
and  ARMAX  models  as  the  most  accurate  models  validates  the  assumption  that  these 
disturbances  to  the  system  entered  the  process  early  on.  The  disturbance  created  by  the 
excitation  of  an  actuator  in  turn  causes  excitation  of  other  actuators,  especially  for 
actuators  that  share  the  same  aluminum  base  in  the  region  above  200  Hz.  Since  the 
disturbance  is  caused  by  another  actuator,  the  choice  of  modeling  the  noise  with  the 
system  dynamics  (poles)  was  an  accurate  one.     The  interaction  or  coupling  among 
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actuators,  particularly  between  base  adjoining  actuators  is  distinct.  Although  the  passive 
damping  mechanism  works  well  in  eliminating  low  frequency  interaction,  high  frequency 
disturbances  are  readily  transmitted  between  actuators.  For  actuators  not  sharing  a 
common  base  node  this  interaction  is  minimal  and  could  possibly  be  ignored.  However, 
for  base  node  sharing  elements  this  problem  is  significant. 

The  lead-lag  compensator  worked  well  for  the  single  actuator  in  the  1Hz  to  100Hz 
range,  providing  over  20dB  damping  of  the  42Hz  disturbance.  However,  above  200  Hz 
the  controller  exacerbates  the  disturbances.  This  problem  is  amplified  in  the  case  where 
all  six  actuators  are  controlled  simultaneously.  The  interaction  between  actuators 
compounds  the  high  frequency  degraded  performance  of  the  controller.  This  causes  it  to 
go  unstable.  The  high  frequency  interaction  between  base  sharing  nodes  must  be 
addressed  before  attempting  to  control  all  six  actuators  simultaneously. 

B.         RECOMMENDATIONS  FOR  FUTURE  WORK 

With  accurate  plant  models  of  the  actuators  the  next  logical  step  is  control  law 
design.  With  knowledge  of  the  system  it  is  possible  to  design  controllers  using  either 
classical  or  modern  techniques.  The  main  concern  in  the  design  is  the  high  frequency 
interaction  in  the  range  of  1kHz  to  1.5  kHz.  The  revelation  of  actuator  interaction  does 
not  necessarily  mean  that  it  is  impossible  to  control  each  actuator  independently.  The 
dynamic  compensator  applied  to  an  individual  actuator  worked  well  below  100Hz  but  did 
not  roll  off  sufficiently.  This  resulted  in  high  frequency  instability  in  control  of  all  six 
actuators.  With  knowledge  of  the  plant  transfer  function  it  is  possible  to  design  a 
dynamic  compensator  capable  of  controlling  each  actuator  independently.  However,  the 
design  must  be  such  that  it  does  not  excite  the  1 .4  kHz  resonance.  On  the  other  hand,  it 
may  be  necessary  to  treat  the  UQP  as  a  multiple  input/multiple  output  (MIMO)  system 
and  a  modern  control  design  is  therefore  required.  Another  possible  remedy  to  the 
interaction  between  base  node  sharing  elements  is  to  change  the  material  of  the 
attachment  node.  The  use  of  aluminum  acts  as  an  excellent  conductor  of  vibration  from 
one  actuator  to  another.  By  using  vibration  absorbing  material  such  as  high  stiffness 
rubber  it  may  be  possible  to  reduce  the  interaction  sufficiently  such  that  each  actuator 
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could  be  controlled  independently  using  single  input/single  output  (SISO)  classical  design 
techniques.  Replacing  the  upper  node  material  could  also  reduce  the  overall  transmission 
of  unwanted  vibration,  not  only  from  interaction  but  external  sources  as  well. 
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APPENDIX  A.  COMPUTER  PROGRAMS 

All  computer  codes  included  in  this  thesis  were  written  by  Naval  Postgraduate 
School  Space  Research  and  Design  Center  unless  otherwise  noted.  The  code  included  in 
this  appendix  is  organized  along  lines  of  the  thesis. 

A.  SYSTEM  IDENTIFICATION  AND  MODEL  CONSTRUCTION 
%  sysidentm  :  MAIN  PROGRAM 

%  Written  by:  George  D.  Beavers  June  1997 

%  This  MATLAB  ".m"  file  reads  in  the  user  requested  actuator  input  output  data 

%  from  the  "4h*.mat"  files.  "*"  represents  the  actuator  that  is  being  excited  by  the 

%  50mV  "pink  noise".  It  calls  the  user  written  functions 

%  "arxmod.m",  "armaxmod.m"  and  "bjmod.m"  to  calculate  models  based  on  the 

%  ARX,  ARMAX  and  Box-Jenkins  parameter  estimation  methods  respectively 

%  The  zeros,  poles  and  transfer  functions  are  then  extracted  from  the  models 

%  returned  by  the  user  written  functions. 

%  Each  "4h*.mat"  file  has  7  rows  of  data 

%         rows  1-6  =  geophone  output  from  the  respective  actuator  (row  1  =actuator  #1) 

%         row  7  is  the  50mv  "pink  noise"  input 

%  Information  on  the  MATLAB  SITB  functions  can  be  found  in  reference[12] 

clear 

dum='y'; 

while  (dum=='y'), 

j=input('Enter  strut  number  »  '); 

disp(['load  4h'  int2str(j)]) 

eval(['load  '  '4h'  int2str(j)  ])  %Load  input  output  data  % 

z2=[trace_y(j,  1:10000)'  trace_y(7, 1:10000)'];  %Pull  out  requested  channel  % 

z2=dtrend(z2);  %Remove  trends  &  adjust  levels% 

z3=[trace_y(j,  1 000 1 :20000)'  trace_y(7, 1 000 1 :20000)'] ; 

z3=dtrend(z3); 

disp(['Choose  Structure  Selection  Method']) 

disp(['l=Automatic'])  %  Use  algorithm  to  determine  % 

disp(['2=Manual'])  %[na,  nb,  nk]  or  input  them      % 

temp=input('»> ');  %  manually  % 


'O 


disp(['Choose  a  Parametric  Estimation  Model']) 
disp(['arx     =1';       'armax  =2';       'bj      =3']) 
chc=input('»>  ') 

%  The  following  subroutine  is  used  to  determine  input  delay  and  number  of  parameters% 
%  and  structure  na,nb  based  on  that  delay  for  an  ARX  based  model  % 
if  temp==  1 

V=arxstruct(z2,z3,struc(2,2, 1 :5)); 
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[nn,Vm]=selstruc(V,'AIC) 
inp=input('enter  delay') 
V=arxstruct(z2,z3,struc(  1:15,1: 15, inp)); 
nnA=selstruct(V,'AIC) 
nnM=selstruct(V,'MDL') 
Figure(l) 

nns=selstruct(V)  %  Generate  parameter  Vs.  Loss  fcn% 

%  plot  % 

%  Determine  which  model  to  construct  then  call  the  function  % 
%  Using  the  automatically  generated  ARX  based  Structure  % 
if  chc==l 

arxmod(j,nns,z2,z3); 
elseif  chc==2 

nc=input('Enter  Order  of  Noise  model  nc  »'); 

itr=input('Enter  Maximum  Number  of  Iterations  »'); 

narm=[nns(l,l)  nns(l,2)  nc  nns(l,3)]; 

armaxmod(j,narm,itr,z2,z3) 


else 


on=input('Enter  Order  of  Noise  model  [nc  nd]  »') 
itr=input('Enter  Maximum  Number  of  Iterations  »') 
nbj=[nns(l,2)  on  nns(l,l)  nns(l,3)] 
bjmod(j,nbj,itr,z2,z3) 


end 


%  Generate  models  based  on  manually  entered  structure  % 
else 

if  chc==l 

nns=input('Enter  Structure  [na  nb  nk]  »'); 
[arxm,arxmft]=arxmodl(j,nns,z2,z3); 
[zepo  1  ,k  1  ]=th2zp(arxm)  %  Extracts  zeros-poles 

[num  1  ,den  1  ]=th2tf(arxm)  %  Extracts  transfer  function 

[zepo  1  f,k  1  f]=th2zp(arxmft)    %  Extracts  zeros-poles  (filtered  data) 
[num  1  f ,den  1  f]=th2tf(arxmft)    %  Extracts  tranfer  function(filtered  data) 
present(arxm)  %  Extracts  parameters 

present(arxmft)  %  Extracts  parameters(filtered  data) 

elseif  chc==2 

narm=input('Enter  Structure  [na  nb  nc  nk]  »'); 
itr=input('Enter  Maximum  Number  of  Iterations  »') 
[armxm,armxmft]=armaxmdf(j,narm,itr,z2,z3); 
[zepo2,k2]=th2zp(armxm)      %  Extracts  zeros-poles 
[num2,den2]=th2tf(armxm)     %  Extracts  tranfer  function 
[zepo2f,k2f]=th2zp(armxmft)  %  Extracts  zeros-poles  (filtered  data) 
[num2f,den2f]=th2tf(armxmft)  %  Extracts  tranfer  function(filtered  data) 
present(armxm)  %  Extracts  parameters 
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else 


end 
end 


present(armxmft)  %  Extracts  parameters(filtered  data) 

nbj=input('Enter  Structure  [nb  nc  nd  nf  nk]  »') 
itr=input('Enter  Maximum  Number  of  Iterations  »') 
[bjm,bjmft]=bjmdf(j,nbj,itr,z2,z3); 


dum=input('Continue  ?» ','s'); 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  arxmod.m  % 

%  ARX  Model  Construction  and  Analysis  % 

%  George  D.  Beavers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  generates  a  model  based  on  the  input-output  data  of  a  system  using  % 
%  the  ARX  parameter  estimation  method.  It  also  generated  a  series  of  graphical  and  % 
%  numerical  results  in  order  to  evaluate  the  suitability  of  the  model  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function  [arxm,arxmft]=arxmod(j,nns,z2,z3) 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  arxm=ARX  based  model  % 

%  arxmft=ARX  based  model  (filtered  data)  % 

%  j=actuator  % 

%  nns=structure([na  nb  nc])  % 

%  z2=input-output  data  for  model  construction        % 
%  Z3=input-ouput  data  for  model  evaluation  % 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

na=nns(l,l); 
NA=int2str(na); 
nb=nns(l,2); 
NB=int2str(nb); 

nk=nns(l,3); 
NK=int2str(nk); 

or  *****************  ARX  MODEL  *************************** 

arxm=arx(z2,nns); 

arxm=sett(arxm,0.  le-3);  %  Set  sample  interval  % 

garxm=th2ff(arxm);  %  Convert  to  frequency  function% 
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%  This  section  loads  separate  data  to  compute  the  actual  actuator  frequency  response% 
%  and  then  plot  it  against  the  model's  frequency  response  % 

Figure(2) 

eval(['load  s',int2str(j),'m'])  % 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

semilogxCX^Yl/b'), 

hold  on 

semilogx(garxm(:,l)/(2*pi),20*logl0(garxm(:,2)),,r'); 

title([,Blue=Actuator#',num2str(j),'Red=ARX  MODEL[\NA,'  ',NB,'  ',NK,T]) 

xlabel('Frequency  (Hz)') 

ylabelCMagnitudeCdB)*) 

grid 

hold  off 

Figure(3) 

eval(['load  s,,int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 

semilogx(X2,Y2,,b') 

hold  on 

semilogx(garxm(:,l)/(2*pi),(garxm(:,4)),'r'); 

title([,Blue=Actuator#,,num2str(j),',  Red=ARX  MODEL[',NA,'  \NB,'  ',NK,T]) 

xlabel('Frequency  (Hz)') 

ylabel(*Phase  (Deg)') 

grid 

hold  off 

%************#*****^j^X  MODEL  (DATA  PR EFILTERED^*********** ******% 

%  The  above  process  is  repeated  using  prefiltered  data.  A  10th  order  Butterworth  filter  % 

%  is  used.  The  passband  is  given  in  terms  of  fractions  of  the  Nyquist  frequency  (0.52)% 

%  This  is  a  low  pass  filter  with  an  upper  bound  of  approximately  1.5kHz 

% 

zf=idfilt(z2, 10,.52);  %Prefilter  model  construction  data% 

arxmft=arx(zf,nns); 
arxmft=sett(arxmft,0.  le-3); 
garxmf t=th2ff (arxmft) ; 

eval(['load  s',int2str(j),'m']) 
Xl=o2ilx; 
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Yl=20.*logl0(abs(o2il)); 

Figure(4) 

semilogx(Xl,Yl,'b') 
hold  on 

semilogx(garxmft(:,l)/(2*pi),20*logl0(ga^xmft(:,2)),,r,); 

title(['Blue=Actuator#,,num2str(j),',  Red=ARX  MODEL[',NA,'  ',NB,'  ',NK,']  *Data 

Prefiltered*']) 

xlabel('Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold  off 

Figure(5) 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 

semilogx(X2,Y2,,b') 

hold  on 

semilogx(garxmft(:,  1  )/(2*pi),(garxmft(:,4)),'r'); 

title(['Blue=Actuator#',num2str(j),',  Red=ARX  MODEL[',NA,'  ',NB,'  \NK,']  *Data 

Prefiltered*']) 

xlabel('Frequency  (Hz)') 

ylabel('Phase  (Deg)') 

grid 

hold  off 

%******************GRAPHICAL  ^j^  NUMERICAL  ANALYSIS***********% 

%  This  section  generates  plots  and  numerical  results  used  in  evaluating  the  model  % 

Figure(6) 

%  Apply  the  input  row  of  the  validation  data  to  the  model  and  then  compare  with       % 

%  the  actual  output  of  the  system  corresponding  to  the  validation  data  and  plot  results% 

[y  h  ,fit]  =compare(z3  ,arxm) ; 

axis([5000,5 100,-1.0,1.0]) 

title(['Plot  of  ARX  MODEL  [',NA,'  ',NB,'  ',NK,']  Output  vs.  Actual  output']) 

text(5020,0.75,'Red=Model  Output,  Blue=  Measured  Output') 

xlabel('Sample') 

ylabel('Output  Magnitude') 

Figure(7) 

[yhf,fitf]=compare(zfv,arxmft);  %Prefiltered  data% 

title(['Plot  of  ARX  MODEL  [',NA,'  ',NB,'  ',NK,']  Output  vs.  Actual  output  *Data 

Prefiltered*']) 
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axis([5000,5100,-1.0,1.0]) 

text(5020,0.75,'Red=Model  Output,  Blue=  Measured  Output') 

xlabel('Sample') 

ylabel('Output  Magnitude') 

%  Calculate  the  residuals  and  plot  the  auto  and  cross  correlation  functions  % 

Figure(8) 

el=resid(z3,arxm);  %  Non-filtered  data 

Figure(9) 

e2=resid(zfv,arxmft);  %  Prefiltered  data 

%  Plot  the  residuals  % 

Figure(lO) 

plot(el,'m'),grid 

title(['Residual  of  ARX  Model  [',NA,'  ',NB,'  ',NK,']']) 

ylabel('Residual') 

xlabel('Sample') 

axis(  [5000,5 100,-0. 1,0.1]) 

Figure(ll) 

plot(e2,'m'),grid 

title(['Residual  of  ARX  Model  [',NA,'  ',NB,'  ',NK,']  *Data  Prefiltered*']) 

axis([5000,5 100,-0.05,0.05]) 

xlabel('Sample') 

ylabel('Residual') 

%  Calculate  Numerical  Results% 

elpos=abs(el); 

e2pos=abs(e2); 

ypos=abs(z3(:,l)); 

ave=mean(elpos)  %  Average  residual 

aveft=mean(e2pos)  %  Average  residual 

outave=mean(ypos)  %  Average  actual  output 

fit  %Mean  Square  Fit  calculated  above 

fitf  %Mean  Square  Fit  calculated  above 

%  Calculate  and  plot  the  actual  output-model  output 

z3out=z3(:,l); 

fori=5000:5100; 

difft(i)=z3out(i)-yhf(i); 
end 

z3out=z3(:,l); 
forj=5000:5100; 

diff(j)=z3outa)-yh(j); 
end 

Figure(14) 
plot(diff,'r') 
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axis([5000,5 100,-0.5,0.5]) 

titletf'Actual  Output  -  ARX  Model  [',NA,'  ',NB,'  ',NK,']  Model']) 

xlabel('Sample') 

Figure(15) 

plot(difft,'r') 

axis([5000,5 100,-0.5,0.5]) 

titletf'Actual  Output  -  ARX  Model  [',NA,'  ',NB,'  ',NK,']  Model  *Data  Prefiltered*']) 

xlabel('Sample') 

avediff=mean(abs(diff))  %Ave.  difference  between  actual  and  model  outputs 

avedft=mean(abs(difft))  %  Ave.  difference  between  actual  and  model  outputs  (filter) 

%  Generate  Zero-Pole  plots 

Figure(12) 

zpplot  1  (th2zp(arxm),3,[],  1 .2) 

titletf'Zero  Pole  Plot  ARX  Model  [',NA,'  ',NB,'  ',NK,']  Model']) 

Figure(13) 

zpplot  1  (th2zp(arxmft)  ,3 ,  [] ,  1 .2) 

titletf'Zero  Pole  Plot  ARX  Model  [',NA,'  ',NB,'  ',NK,']  Model  *Data  Prefiltered*']) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  armaxmod.m  % 

%  ARMAX  Model  Construction  and  Analysis  % 

%  George  D.  Beavers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  generates  a  model  based  on  the  input-output  data  of  a  system  using  % 
%  the  ARMAX  parameter  estimation  method.  It  also  generated  a  series  of  graphical  % 
%  and  numerical  results  in  order  to  evaluate  the  suitability  of  the  model  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function  [arxm,arxmft]=armaxmod(j,nns,itr,z2,z3) 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  arxm=ARMAX  based  model  % 

%  arxmft=ARMAX  based  model  (filtered  data)        % 
%  j=actuator  % 

%  nns=structure([na  nb  nc])  % 

%  z2=input-output  data  for  model  construction         % 
%  Z3=input-ouput  data  for  model  evaluation  % 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

na=nns(l,l); 
NA=int2str(na); 
nb=nns(l,2); 
NB=int2str(nb); 
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nc=nns(l,3); 
NC=int2str(nc); 
nk=nns(l,4); 
NK=int2str(nk); 


%*##**#**#****#***^j^l^X  MODEL  **********************************% 
%  This  section  loads  separate  data  to  compute  the  actual  actuator  frequency  response% 
%  and  then  plot  it  against  the  model's  frequency  response  % 

arxm=armax(z2,nns,itr); 

arxm=sett(arxm,0.1e-3);  %  Set  sampling  interval  % 

garxm=th2ff(arxm);  %  Convert  to  frequency  function% 


eval(['load  s',int2str(j),'m'])  %  Load  actual  actuator  data  (magnitude)% 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

Figure(2) 

semilogx(Xl,Yl,*b'), 

hold  on 

semilogx(garxm(:J)/(2*pi),20*loglO(garxm(:,2)),'r'); 

title(['Blue=Actuator#,,num2str(j),*Red=ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  ',NK,T]) 

xlabel('Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold  off 

eval(['load  s',int2str(j),'p'])  %  Load  actual  actuator  data  (phase)% 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 

Figure(3) 

semilogx(X2,Y2,'b*) 

hold  on 

semilogx(garxm(:,l)/(2*pi),(garxm(:,4)),'r'); 

title([,Blue=Actuator#',num2str(j),*Red=ARMAX  MODEL[',NA,'  \NB,'  \NC/  ',NK,']']) 

xlabel('Frequency  (Hz)') 

ylabel('Phase  (Deg)') 

grid 

hold  off 

%*************^j^/|^X  MODEL  (DATA  pR£FILTERED********************% 
%  The  above  process  is  repeated  using  prefiltered  data.  A  10th  order  Butterworth  filter  % 
%  is  used.  The  passband  is  given  in  terms  of  fractions  of  the  Nyquist  frequency  (0.52)% 
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%  This  is  a  low  pass  filter  with  an  upper  bound  of  approximately  1 .5kHz 

% 

zf=idfilt(z2,10,.52);  %  Prefilter  model  construction  data  % 

arxmft=armax(zf,nns,itr); 

arxmft=sett(arxmft,0.  le-3); 

garxmf t=th2ff (arxmft) ; 

eval(['load  s',int2str(j),'m']) 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

Figure(4) 

semilogx(X  1 ,  Y 1  ,'m') 

hold  on 

semilogx(garxmft(:,l)/(2*pi),20*logl0(garxmft(:,2)),,g'); 

titie(['Blue=Actuator#',num2str(j),'Red=ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  'J4K,'] 

*Data  Prefiltered*']) 

xlabel( 'Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold  off 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*un  wrap(angle(o2i  1  ))/pi; 

Figure(5) 

semilogx(X2,Y2,'b') 

hold  on 

semilogx(garxmft(:,l)/(2*pi),(garxmft(:,4)),'r'); 

title(['Blue=Actuator#',num2str(j),'Red=ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  ',NK,'] 

*Data  Prefiltered*']) 

xlabel('Frequency  (Hz)') 

ylabel('Phase  (Deg)') 

grid 

hold  off 

%******************qraPHICAL  fr^jy  NUMERICAL  ANALYSIS***********% 

%  This  section  generates  plots  and  numerical  results  used  in  evaluating  the  model         % 

%  Apply  the  input  row  of  the  validation  data  to  the  model  and  then  compare  with  % 

%  the  actual  output  of  the  system  corresponding  to  the  validation  data  and  plot  results 

% 

Figure(6) 

[yh,fit]=compare  1  (z3,arxm); 

axis([5000,5 100,-1.0,1-0]) 
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title(['Plot  of  ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  ',NK,']  Output  vs.  Actual  output']) 

text(5020,0.75,'Red:Model  Output,  Green:  Measured  Output') 

grid 

Figure(7) 

[yhf,fitf]=compare  1  (zfv,arxmft); 

titletf'Plot  of  ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  \NK,']  Output  vs.  Actual  output 

*Data  Prefiltered*']) 

axis([5000,5100,-1.0,1.0]) 

xlabel('Sample') 

grid 

%  Calculate  residuals  and  generate  and  plot  auto  &  cross  correlation  functions  % 

Figure(8) 

el=residl(z3,arxm); 

Figure(9) 

e2=resid  1  (zfv,arxmft); 

%  Plot  Residuals    % 

Figure(lO) 

plot(el,'m'),grid 

title(['Residual  of  ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  ',NK,']']) 

ylabel('Residual') 

axis([5000,5100,-0.1,0.1]) 

Figure(ll) 

plot(e2,'m'),grid 

title(['Residual  of  ARMAX  MODEL[',NA,*  ',NB,'  ',NC,'  ',NK,']  *Data  Prefiltered*']) 

axis([5000,5 100,-0.05,0.05]) 

xlabel('Sample') 

ylabel('Residual') 

%  Calculate  Numerical  Results% 

elpos=abs(el); 

e2pos=abs(e2); 

ypos=abs(z3(:,l)); 

ave=mean(elpos)  %  Average  residual 

aveft=mean(e2pos)  %  Average  residual  (Prefiltered) 

outave=mean(ypos)  %  Average  actual  output 

fit  %  Mean  Square  Fit  calculated  above 

fitf  %  Mean  Square  Fit  (Prefiltered) 

%  Calculate  actual  output  minus  model  output  an  plot  results  % 

z3out=z3(:,l); 
fori=5000:5100; 
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difft(i)=z3out(i)-yhf(i); 
end 

z3out=z3(:,l); 
forj=5000:5100; 

diffO)=z3out(j)-yh(j); 
end 

Figure(14) 
plot(diff,'r') 

axis([50O0,5 100,-0.5,0.5]) 

title(['Actual  Output  -  ARMAX[',NA,'  ',NB,'  ',NC,'  ',NK,*]  Model  Output']) 
grid 

Figure(15) 
plotCdifft/r1) 

axis([5000,5 100,-0.5,0.5]) 

title(['Actual  Output  -  ARMAX[',NA,'  ',NB,'  ',NC,'  ',NK,']Model  Output  *Data 
Filtered*']) 
grid 

avediff=mean(abs(diff))  %  Ave.  output  difference 

avedft=mean(abs(difft))  %  Ave.  output  difference  (Prefiltered) 

%  Generate  Zero-Pole  Plots 

Figure(12) 

zpplot(th2zp(arxm),3,[],1.2) 

title(['Zero  Pole  Plot  ARMAX  MODEL[',NA,'  ',NB,'  ',NC,'  ',NK,']']) 

Figure(13) 

zpplot(th2zp(arxmft),3,[],1.2) 

title(['Zero  Pole  Plot  ARMAX  [',NA,'  ',NB,'  ',NC,'  ',NK,']  Model  *Data  Prefiltered*']) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  bjmod.m  % 

%  Box-Jenkins  Model  Construction  and  Analysis  % 

%  George  D.  Beavers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  generates  a  model  based  on  the  input-output  data  of  a  system  using  % 
%  the  Box-Jenkins  parameter  estimation  method.  It  also  generated  a  series  of  graphical% 
%  and  numerical  results  in  order  to  evaluate  the  suitability  of  the  model  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function  [bjm,bjmft]=bjmod(j,nns,itr,z2,z3) 
%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  bjm=Box-Jenkins  based  model  % 

%  bjmft=Box-Jenkins  based  model  (filtered  data)     % 
%  j=actuator  % 

%  nns=structure([na  nb  nc])  % 

%  itr=number  of  iterations  % 
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%  z2=input-output  data  for  model  construction        % 

%  Z3=input-ouput  data  for  model  evaluation  % 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

nb=nns(l,l); 

NB=int2str(nb); 

nc=nns(l,2); 

NC=int2str(nc); 

nd=nns(l,3); 

ND=int2str(nd); 

nf=nns(l,4); 

NF=int2str(nf); 

nk=nns(l,5); 

NK=int2str(nk); 

%*****************gQX_j^jsJKINS  ^q£)£l*************************  ******% 
%  This  section  loads  separate  data  to  compute  the  actual  actuator  frequency  response  % 
%  and  then  plot  it  against  the  model's  frequency  response  % 

arxm=bj(z2,nns,itr);  %Build  the  model  % 

arxm=sett(arxm,0.1e-3);  %  Set  the  sampling  interval  % 

garxm=th2ff(arxm);  %  Convert  to  frequency  function        % 

eval(['load  s',int2str(j),'m'])  %Load  data  for  actuator  frequency  function% 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

Figure(2) 

semilogx(X  1 ,  Y 1  ,'b'), 

hold  on 

semilogx(garxm( : ,  1  )/(2*pi),20*log  1 0(garxm(:  ,2)),'r'); 

title(['Blue=Actuator#',num2str(j),'Red=Box-Jenkins  MODEL[',NB,'  ',NC,'  ',ND,'  ',NF,' 

',NK,T]) 

xlabel('Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold  off 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 

Figure(3) 

semilogx(X2,Y2,'b') 

hold  on 

semilogx(garxm(:,  1  )/(2*pi),(garxm(:,4)),'r'); 

title(['Blue=Actuator#',num2str(j),,Red=Box-Jenkins  MODEL[',NB,'  ',NC,'  ',ND,'  ',NF,' 

\NK,T]) 
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xlabel('Frequency  (Hz)') 
ylabel('Phase  (Deg)') 
grid 
hold  off 

%*************BOx-JENKINS  MODEL  (DATA  PREFTLTERED**************% 
%  The  above  process  is  repeated  using  prefiltered  data.  A  10th  order  Butterworth  filter  % 
%  is  used.  The  passband  is  given  in  terms  of  fractions  of  the  Nyquist  frequency  (0.52)% 
%  This  is  a  low  pass  filter  with  an  upper  bound  of  approximately  1 .5kHz  % 

zf=idfilt(z2,10,.52);  %  Prefilter  model  construction  data  % 

arxmft=bj(zf,nns,itr);  %  Construct  model  % 

arxmft=sett(arxmft,0.1e-3);  %  Set  sampling  interval  % 

garxmft=th2ff(arxmft);  %  Convert  to  frequency  function       % 

eval(['load  s',int2str(j),'m'])  %  Load  data  for  actuator  frequency  response 

Xl=o2ilx; 

Yl=20.*logl0(abs(o2il)); 

Figure(4) 

semilogx^UYl/b') 

hold  on 

semilogx(garxmft(:,  1  )/(2*pi),20*log  10(garxmft(:,2)),'r'); 

title([,Blue=Actuator#*,num2str(j),'Red=Box-Jenkins  MODEL[',NB,'  ',NC,'  ',ND,'  ',NF,' 

*,NK,*]*Data  Prefiltered*']) 

xlabel('Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

grid 

hold  off 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 

Figure(5) 

semilogx(X2,Y2,,m') 

hold  on 

semilogx(garxmft(:,l)/(2*pi),(garxmft(:,4)),'g'); 

title([,Blue=Actuator#',num2str(j),,Red=Box- Jenkins  MODEL[',NB,"  \NC,'  ',ND,'  ',NF,' 

',NK,']*Data  Prefiltered*']) 

xlabelfFrequency  (Hz)') 

ylabel('Phase  (Deg)') 

grid 

hold  off 
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%******************GRAPHICAL  AND  NUMERICAL  ANALYSIS***********% 
%  This  section  generates  plots  and  numerical  results  used  in  evaluating  the  model  % 
%  Comparison  of  model  output  to  actual  output  % 

%  Apply  the  input  row  of  the  validation  data  to  the  model  and  then  compare  with  % 

%  the  actual  output  of  the  system  corresponding  to  the  validation  data  and  plot  results 
% 

Figure(6) 

[yh,fit]=compare(z3,arxm); 

axis([5000,5100,-1.0,1.0]) 

title([*Plot  of  Box-Jenkins  [',NB,'  ',NC,'  \ND,'  \NF,'  ',NK,']  Output  vs.  Actual  output']) 

text(5020,0.75,'Red:Model  Output,  Blue:  Measured  Output') 

xlabel('Sample') 

Figure(7) 

[yhf ,f itf]=compare(z3  ,arxmft) ; 

title(['Plot  of  Box-Jenkins  MODEL  [',NB,'  ',NC,'  ',ND,'  ',NF,'  ',NK,']  Output  vs.  Actual 

output  *Data  Prefiltered*']) 

text(5020,0.75,'Red:Model  Output,  Blue:  Measured  Output') 

axis([5000,5100,-1.0,1.0]) 

xlabel('Sample') 

%  Calculate  residuals  and  auto  and  cross  correlation  functions  and  plot  results 

Figure(8) 

e  1  =resid  1  (zfv,arxm); 

Figure(9) 

e2=residl(zfv,arxmft); 

Figure(lO) 

plot(el,'m'),grid 

title( ['Residual  of  Box-Jenkins  Model  [',NB,'  ',NC,'  ',ND,'  ',NF,*  ',NK,']']) 

ylabel('Residual') 

xlabel('Sample') 

axis([5000,5 100,-0. 1,0.1]) 

Figure(ll) 

plot(e2,'m'),grid 

title(['Residual  of  Box-Jenkins  Model  [',NB,'  ',NC,'  ',ND,'  ',NF,'  ',NK,']  *Data 

Prefiltered*']) 

axis([5000,5 100,-0.05,0.05]) 

xlabel('Sample') 

ylabel('Residual') 
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%  Numerical  results 

elpos=abs(el); 

e2pos=abs(e2); 

ypos=abs(z3(:,l)); 

ave=mean(elpos)  %  Average  residual 

aveft=mean(e2pos)  %  Average  residual  (Prefiltered  Data) 

outave=mean(ypos)  %  Average  actuator  output 

fit  %  Mean  Square  Fit  calculated  above 

fitf  %  Mean  Square  Fit  (Prefiltered) 

%  Calculate  and  plot  actual  output  minus  model  output 

z3out=z3(:,l); 

for  1=5000:5100; 

difft(i)=z3out(i)-yhf(i); 
end 

z3out=z3(:,l); 
forj=5000:5100; 

diff(j)=z3out(j)-yhG); 
end 

Figure(12) 
plot(diff,'r') 

axis([5000,5 100,-0.5,0.5]) 
title('Actual  Output  -  Model  Output') 
Figure(13) 
plot(difft,'r') 

axis([5000,5 100,-0.5,0.5]) 

title('Actual  Output  -  Model  Output  *Data  Filtered*') 
avediff=mean(abs(diff))  %  Average  output  difference 

avedft=mean(abs(difft))  %  Average  output  difference(Prefiltered) 

%  Generate  Zero-Pole  Plot 

Figure(14) 

zpplotl(th2zp(arxm),3,[],1.2) 

title(['Zero  Pole  Plot  Box-Jenkins  Model  [',NB,'  ',NC,'  ',ND,'  ',NF,'  ',NK,'] ']) 

Figure(15) 

zpplot  1  (th2zp(arxmft)  ,3 ,  [] ,  1 .2) 

title(['Zero  Pole  Plot  Box-Jenkins  Model  [',NB,'  ',NC,'  ',ND,'  ',NF,'  ',NK,']  *Data 

Prefiltered*']) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  plotarx.m  % 

%  Compare  Three  ARX  Model  Structures  % 

%  George  D.  Beavers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%  This  program  generates  the  frequency  response  plots  of  the  three  best  ARX  model     % 
%  structures  for  evaluation  and  structure  selection  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


j=input('enter  strut»'); 

eval(['load  4h',int2str(j)]) 
eval(['load  s',int2str(j),'m']) 

z2=[trace_y(j,  1 :  10000)'  trace_y(7, 1 :  10000)']; 

z2=dtrend(z2); 

z3=[trace_y (j,  1 000 1 :20000)'  trace_y (7, 1 000 1 :20000)'] ; 

z3=dtrend(z3); 

nnl=[7  6  1]  %  Top  three  structures 

nn2=[15  15  1] 
nn3=[12  12  1] 

Xl=o2ilx; 
Yl=20.*logl0(abs(o2il)); 


zf=idfilt(z2,10,.52); 

%  Construct  the  models  % 

fml=arx(zf,nnl); 

fml=sett(fml,0.1e-3); 

gfml=th2ff(fml); 

fm2=arx(zf,nn2); 

fm2=sett(fm2,0.1e-3); 

gfm2=th2ff(fm2); 


fm3=arx(zf,nn3); 

fm3=sett(fm3,0.1e-3); 

gfm3=th2ff(fm3); 

%  Generate  comparison  plots 

Figure(  1 ) 

semilogx(Xl,Yl,'b'), 

title(['Strut  #',int2str(j),'  Magnitude']) 

xlabel('Frequency  (Hz)') 

ylabel('Magnitude(dB)') 

hold  on 
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semilogx(gfml(:,l)/(2*pi),20*logl0(gfml( 
semilogx(gfm2(:,  1  )/(2*pi),20*log  1 0(gfm2( 
semilogx(gfm3(:,l)/(2*pi),20*loglO(gfm3( 
hold  off 

Figure(2); 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*un  wrap(angle(o2i  1  ))/pi; 

semilogx(X2,Y2,,b'), 

title(['Strut  #',int2str(j),'  Phase']) 

xlabel('Frequency  (Hz)'); 

ylabel('Deg'); 

hold  on 


,2)),'r'); 

,2)),'m'); 

,2)),'w');,gnd; 


semilogx(gfml  (:,  1  )/(2*pi),(gfml  ( 
semilogx(gfm2(:,  1  )/(2*pi),(gfm2( 
semilogx(gfm3(:,  1  )/(2*pi),(gfm3( 


,4)),'r'); 

,4)),'m'); 

,4)),'w');,grid; 


hold  off 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  plotamb.m  % 

%  Model  Comparison  % 

%  George  D.  Beavers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  generates  the  frequency  response  plots  of  the  three  best  models  based  % 
%  on  the  ARX  structure^  6  1].    ARX,  ARMAX  and  BOX- JENKINS  models  are  % 

%  compared  to  select  the  best  model  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
j=input('enter  strut»'); 

eval(['load  4h',int2str(j)]) 
eval(['load  s',int2str(j),'m']) 

z2=[trace_yG,l :  10000)'  trace_y(7,l :  10000)']; 

z2=dtrend(z2); 

z3=[trace_y (j ,  1 000 1 :20000)'  trace_y (7, 1 000 1 :20000)'] ; 

z3=dtrend(z3); 

%  Define  model  structure  % 
nnl=[7  6  2  1] 
nn2=[7  6  1] 
nn3=[62  2  7  1] 

Xl=o2ilx; 
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Yl=20.*logl0(abs(o2il)); 

%  Generate  models  (non-filtered)  % 

%m  1 =armax(z2,nn  1 ,60) ; 

%ml=sett(ml,0.1e-3); 

%gml=th2ff(ml); 

%m2=arx(z2,nn2); 

%m2=sett(m2,0.  le-3); 

%gm2=th2ff(m2); 

%m3=bj(z2,nn3,60); 

%m3=sett(m3,0.1e-3); 

%gm3=th2ff(m3); 


zf=idfilt(z2,10,.52); 

%  Generate  Models  (Prefiltered) 
fm  1  =armax(zf  ,nn  1 ,60) ; 
fml=sett(fml,0.1e-3); 
gfml=th2ff(fml); 

fm2=arx(zf,nn2); 

fm2=sett(fm2,0.1e-3); 

gfm2=th2ff(fm2); 


fm3=bj(zf,nn3,60); 

fm3=sett(fm3,0.1e-3); 

gfm3=th2ff(fm3); 

%  Generate  plots 

Figure(l) 

semilogx(Xl,Yl/b') 

title(['Strut  #',int2str(j),'  Magnitude']) 

xlabel('Frequency  (Hz)') 

y  label('Magnitude(dB )') 

hold  on 

semilogx(gfml(:,l)/(2*pi),20*logl0(gfml( 

semilogx(gfm2(:,l)/(2*pi),20*logl0(gfm2( 

semilogx(gfm3(:,l)/(2*pi),20*logl0(gfm3( 

hold  off 

Figure(2); 

eval(['load  s',int2str(j),'p']) 

X2=o2ilx; 

Y2=  1 80*unwrap(angle(o2i  1  ))/pi; 


,2)),'r'); 

,2)),'m'); 

,2)),'w');,grid; 
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semilogx(X2,Y2,'b') 
title(['Strut  #',int2str(j),'  Phase']) 
xlabel('Frequency  (Hz)'); 
ylabel('Deg'); 
hold  on 


semilogx(gfml( 
semilogx(gfm2( 
semilogx(gfm3( 


,l)/(2*pi),(gfml( 
,l)/(2*pi),(gfm2( 
,l)/(2*pi),(gfm3( 


,4))/r'); 

,4)),'m'); 

,4)),'w');,grid; 


hold  off 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  compid.m  % 

%  Interaction  Comparison  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  generates  the  frequency  response  plots  of  the  actuators  when  actuator  % 
%  number  one  is  excited  with  50mV  "pink  noise".  These  plots  use  the  ARX[7  6  1]        % 
%  model  frequency  repines  of  the  respective  actuator.  This  is  compared  with  the  % 

%  modeled  response  of  actuator  number  one  These  plots  are  used  to  check  for  % 

%  interaction  between  actuators  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function  []=compid(j) 
%  j=  actuator  number  % 
int2str(j) 

disp(['load  4h'  int2str(j)]) 
eval(['load  '  '4h'  int2str(j)  '.mat']) 

z  1 =[trace_y ( 1,1:1 0000)'  trace_y (7, 1 : 1 0000)'] ; 

zl=dtrend(zl); 

zla=[trace_y(l,  10001:20000)'  trace_y(7, 10001:20000)']; 

zla=dtrend(zla); 

thl=arx(zl,[7  6  1]); 

thl=sett(thl,0.1e-3); 

gthl=th2ff(thl); 


for  i=2:6 


z2=[trace_y (i,  1 : 1 0000)'  trace_y (7, 1 : 1 0000)'] ; 

z2=dtrend(z2); 

z7=[trace_y(i,  10001:20000)'  trace_y(7, 10001 :20000)']; 
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z7=dtrend(z7); 

th=arx(z2,[7  6  1]); 
th=sett(th,0.1e-3); 
gth=th2ff(th); 


Figure(i-l); 

semilogx(gthl(:,l)/(2*pi),20*logl0(gthl(:,2)),,b'); 
hold  on 

semilogx(gth(:,l)/(2*pi),20*logl0(gth(:,2)),y); 
title('Magnitude:  Input=50mV  Pink  Noise') 
xlabel('Frequency  (Hz)') 
ylabel('dB') 

text(  10,45, ['Excite  Actuator  #'  num2str(j)  'Blue;  Sense  Actuator  ', 
num2str(i),'=Red']) 
grid 
hold  off 

end 

for  i=2:6 

z2=[trace_y(i,  1 :  10000)'  trace_y(7, 1 :  10000)']; 

z2=dtrend(z2); 

z7=[trace_y(i,  10001 :20000)'  trace_y(7, 10001 :20000)']; 

z7=dtrend(z7); 

th=arx(z2,[7  6  1]); 
th=sett(th,0.1e-3); 
gth=th2ff(th); 

Figure(i+4) 

semilogx(gthl(:,l)/(2*pi),(gthl(:,4));b'); 
hold  on 

semilogx(gth(:,l)/(2*pi),(gth(:,4)),'r'); 
title('Phase:  Input=50mv  Pink  Noise') 

text(10,-550,['Excite  Actuator  #'  num2str(j)  -Blue;  Sense  Actuator  ',num2str(i), 
'=Red']) 

xlabel('Frequency  (Hz)') 
ylabel('Degrees') 
grid 
hold  off 
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end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  figl.m,  fijj2.ni,  fig3.m,  flg4.m  % 

%  PSD  of  Applied  Controllers  % 

%  June  1997  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  This  program  generates  the  power  spectral  density  plots  resulting  from  applying         % 

%  a  lead-lag  compensator  to  an  individual  actuator  and  the  all  six  actuators  to  evaluate  % 

%  the  effects  of  actuator  interaction  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%figl.m% 

%  Plots  over  range  of  0-5kHz% 

clear 

load  si  %Open  loop  data  with  42  Hz  disturbance  on% 

Figure(l) 

clg 

psd(trace_y(2, 1 :20000),  1024*  1 ,  10000,[]); 

%[popen,fopen]=psd(trace_y(2,l  :20000),  1024*  12, 10000,[]); 

hold  on 

load  s2  %Compensator  applied  to  actuator  #2  ( 

psd(trace_y(2, 1 :20000),  1024*  1 , 1 0000,kaiser(5 1 2,5),256); 

%[ps2,fs2]=psd(trace_y(2,l:20000),1024*12,10000,[]); 

load  s3  %Compensator  applied  to  all  six  actuators 

psd(trace_y(2, 1 :20000),  1 024*  1 ,  10000,[]); 

%  [ps3,fs3]=psd(trace_y(2, 1 :20000),  1024*  1 2, 1 0000,  []); 

hold  off 

%semilogx((fopenk),20*log  1 0(popenk),'r',(fs2k),20*log  1 0(fs2k),'y'  ,(fs3k),20*log  1 0(fs3k) 

,'b1) 

Title('Blue-Open  Loop;  Magenta-Strut  #2  Closed  Loop;  Red-All  6  Struts  Closed  Loop' ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig2.m% 

%  Plots  same  Figure  as  figl.m  over  range  of  0-200  Hz% 

clear 

load  si  %Open  loop-42  Hz  Disturbance  on% 

Figure(2) 

clg 

psdb(trace_y(2, 1:20000),  1024*8, 10000,[]); 

%[popen,fopen]=psd(trace_y(2, 1 :20000),  1 024*  1 2, 1 0000,[]); 

hold  on 
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load  s2  %Compensator  applied  to  actuator  #2% 

psd(trace_y(2, 1 :20000),  1024*8, 1 0000,[]); 

%[ps2,fs2]=psd(trace_y(2,l:20000),1024*12,10000,[]); 

load  s3  %  Compensator  applied  to  actuator  #2% 

psd(trace_y(2,l:20000),1024*8,10000,[]); 

%[ps3,fs3]=psd(trace_y(2,l:20000),1024*12,10000,[]); 

hold  off 

%semilogx((fopenk),20*logl0(popenk),,r',(fs2k),20*logl0(fs2k),*y',(fs3k),20*logl0(fs3k) 

,'b') 

Title('Blue-Open  Loop;  Magenta-Strut  #2  Closed  Loop;  Red-All  6  Struts  Closed  Loop' ) 

%Title('Open  Loop-42Hz  Disturbance' ) 

axis([0  200  -70  30]) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig3.m% 

%plots  PSD  sensed  by  3  geophones;  actuator  #2,  #3  and  #6  over  range  of  0-200Hz% 

clear 

load  s2  %  Compensator  applied  to  actuator  #2  % 

Figure(3) 

clg 

psd(trace_y(2, 1 :20000),  1024*  1 ,10000,[]); 

%[popen,fopen]=psd(trace_y(2,l:20000),1024*12,10000,[]); 

hold  on 

psd(trace_y(3, 1 :20000),  1 024*  1 ,  10000,[]); 

%[ps2,fs2]=psd(trace_y(3,l:20000),1024*12,10000,[]); 

psd(trace_y(6, 1 :20000),  1024*  1 ,  10000,[]); 

%[ps3,fs3]=psd(trace_y(6,l:20000),1024*12,10000,[]); 

hold  off 

%semilogx((fopenk),20*logl0(popenk),'r,,(fs2k),20*logl0(fs2k),y,(fs3k),20*logl0(fs3k) 

/b') 

Title('Blue-#2  Geophone;  Magenta-#3  Geophone;  Red-#  6  Geophone' ) 

%axis([0  200  -70  30]) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fig3.m% 

%plots  PSD  sensed  by  3  geophones;  actuator  #2,  #3  and  #6  over  range  of  0-5kHz% 

clear 

load  s2  %Compensator  applied  to  actuator  #2  % 

Figure(4) 

clg 

psd(trace_y(2, 1 :20000),  1 024*8, 1 0000,[]); 

%[popen,fopen]=psd(trace_y(2,l:20000),1024*12,10000,[]); 

hold  on 
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psd(trace_y(3,l:20000),l  024*8, 10000,[]); 

%[ps2,fs2]=psd(trace_y(3, 1 :20000),  1024*  12, 10000,[]); 

psd(trace_y  (6, 1 :20000),  1 024*8, 1 0000,[]); 

%[ps3,fs3]=psd(trace_y(6, 1:20000),  1024*  12, 10000,[]); 

hold  off 

%semilogx((fopenk),20*loglO(popenk),'r',(fs2k),20*loglO(fs2k),,y',(fs3k),20*loglO(fs3k) 

■V) 

Title('Blue-#2  Geophone;  Magenta-#3  Geophone;  Red-#  6  Geophone' ) 
axis([0  200  -70  30]) 


ill 
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APPENDIX  B.  MODEL  STRUCUTURE  SELECTION  PLOTS 

A.  Structure  Acceptance  Plots  for  Actuator  Two 

The  following  plots  were  utilized  in  selecting  the  ARX  structure  [na=7  nb=6  nk=l] 
for  actuator  number  two. 
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Plot  of  ARX  MODEL  [7  6  1]  Output  vs  Actual  output  'Data  Prefiltered' 
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B.  Model  Structure  Acceptance  for  Actuator  Three 

The  following  plots  were  utilized  in  selecting  the  ARX  structure  [na=7  nb=7  nk=  1  ] 
for  actuator  number  three. 
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Plot  of  ARX  MODEL  [7  6  1]  Output  vs.  Actual  output 
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C.  Model  Structure  Acceptance  for  Actuator  Four 

The  following  plots  were  utilized  in  selecting  the  ARX  structure  [na=7  nb=7  nk=l] 
for  actuator  number  four. 


RETURN  TO  COMMAND  SCREEN  TO  SELECT  ft  OF  PARAMETERS  TO  BE  ESTIMATED 


15  20 

#  of  par's 


Figure  C  1 
Parameter  Number  Vs.  Loss  Function  For  nk=\ 


Zero  Pole  Plot  ARX  Model  [7  6  1]  Model  'Data  Prefiltered* 
1.5 


-0.5 


Figure  C.2 
Zero-Pole  Plot 


121 


60 


40 


Blue=Actuator#4,  Red=ARX  MODEL[7  6  1]  *Data  Prefiltered* 


0Q 
T5 

20 

CD 

-J 

c 

O) 

oj 

0 

-20 


-40 


^* 

<■- 

\ 

\ 

i 

/ 

1 

i 

10L 


10  10 

Frequency  (Hz) 

(a)  Magnitude 


10" 


1<T 


400 


300 


200 


g    100 
to 

i?     o 

CL 


■100 


-200 


Blue=Actuator#4,  Red=ARX  MODEL[7  6  1]  *Data  Prefiltered' 


-300 


\ 

\ 

\ 
\ 

s 

/ 

/ 

"""■^5: 

^^ 

*j, 

I 

u 

>l 

10l 


10  10 

Frequency  (Hz) 

(b)  Phase 


10" 


10* 


Figure  C.3 
Frequency  Response  Comparison 


122 
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D.  Model  Structure  Acceptance  for  Actuator  Five 

The  following  plots  were  utilized  in  selecting  the  ARX  structure  [na=7  nb=7  nk=l] 
for  actuator  number  five. 
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Figure  D.7 
Actual  Output  and  Model  Output  Difference 
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E.  Model  Structure  Acceptance  for  Actuator  Six 

The  following  plots  were  utilized  in  selecting  the  ARX  structure  [na=7  nb=7  nk=l] 
for  actuator  number  six. 


RETURN  TO  COMMAND  SCREEN  TO  SELECT  #  OF  PARAMETERS  TO  BE  ESTIMATED 
0.025 


0.02 


0.015 


0.01 


0.005  - 


10  15  20 

#  of  par's 


25 


Figure  E.  1 
Parameter  Number  Vs.  Loss  Function  For  nk=\ 


Zero  Pole  Plot  ARX  Model  [7  6  1]  Model  *Data  Prefirtered* 
1.5 


-1.5         -1         -0.5  0  0.5  1  1.5 


Figure  E.2 
Zero-Pole  Plot 


129 


50 


40 


30 


20 


Blue=Actuator#6,  Red=ARX  MODEL[7  6  1]  *Data  Prefiltered' 


QQ 

3 


TO 


■10 


-20 


t    FV 

1 

I 

/^' 

\ 

1 

, 

/ 

/'"/ 

.r 

— 

J 

10u 


10  10 

Frequency  (Hz) 

(a)  Magnitude 


10' 


10H 


200 


100 


o>       0 
Q 

CD 

w 

TO 

£-100 


Blue=Actuaior#6,  Red=ARX  MODEL[7  6  1]  *Data  Prefiltered' 


-200 


-300 


A 

v 

/ 

V 

^ 

** 

*»», 

• 

'  \ 

\ 

\i 

h5 

i 

10^ 


10  10  10 

Frequency  (Hz) 

(b)  Phase 

Figure  E.3 
Frequency  Response  Comparison 


10H 


130 


PlotofARXMODEL[76  1]  Output  vs.  Actual  output 'Data  PrefilterecT 
1 


5000 


0.5 


-0.5 


5020  5040  5060 

Sample 


5080 


5100 


Figure  E.4 
Model  Output  Vs.  Actual  Output 

Correlation  function  of  residuals.  Output  #  1 


1                         -~l~ 

I T 

i               i 

0  5  10  15  20  25 

Cross  corr.  function  between  input  1  and  residuals  from  output  1 


Figure  E.5 
Auto  and  Cross  Correlation  Functions 


131 


0.05 
0.04 


0.03 


0.02 

_   0.01 

ro 

1      o 

ID 

a: 


-0.01 
-0.02 
-0.03 
-0.04 


-0  05 


Residual  of  ARX  Model  [7  6  1]  *Data  Prefiltered* 


5000  5020  5040  5060 

Sample 

Figure  E.  6 
Residual  Vs  Sample 


5080 


5100 


Actual  Output  -  ARX  Model  [7  6  1]  Model  *Data  Prefiltered* 


5000  5020  5040  5060 

Sample 


5080 


5100 


Figure  E.7 
Actual  Output  and  Model  Output  Difference 


132 


APPENDIX  C.  MODEL  VALD3ATION  PLOTS 

A.  Application  of  ARMAX  Method  to  Actuator  Number  One 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX  and 
Box-Jenkins  parameter  estimation  methods  in  modeling  actuator  number  one.  The 
structure  [na=7  nb=6  nc=2  nk=l]  was  used 


Blue=Actuator#1Red=ARMAX  MODEL[7  6  2  1]  'Data  Prefiltered* 


60 


40 


B  20 
S 

■D 

B 

c 
(0      0 


-20 


-40 

10u  10'  10'  10'  10* 

Frequency  (Hz) 

(a)  Magnitude 

Blue=Actuator#1  Red=ARMAX  MODEL[7  6  2  1]  *Data  Prefiltered* 
200 


A*\ 

/'    ' 

y 

r> 

\ 
\ 

/ 

/ 

■ 

■ 

1 

/ 

0) 

Q 
^T-200 

(/> 
(D 

JZ 

a. 


^100 


-600 


* 

V 

Jr 

^^ 

iJr- 

) 

\ 

10 


10' 


10' 
Frequency  (Hz) 

(b)  Phase 
Figure  A.  1 


10° 


10' 


133 


Frequency  Response  Comparison 
Plot  of  ARMAX  MODEL[7  6  2  1]  Output  vs.  Actual  output  *Data  Prefiltered' 
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Zero  Pole  Plot  ARMAX  [7  6  2  1]  Model  *Data  Prefiltered* 
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B.  Application  of  ARMAX  Method  to  Actuator  Number  Two 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX 
parameter  estimation  method  in  modeling  actuator  number  two.  The  structure  [na=7  nb=6 
nc=2  nk=l]  was  used  Two  additional  plots  from  the  Box- Jenkins  [nb=6  nc-2  nc=2  na=l 
nk=\]  model  are  included  for  comparison. 
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Figure  B.  1 
Frequency  Response  Comparison 
PlotofARMAXMODEL[7  6  2  1]  Output  vs.  Actual  output 'Data  Prefiltered' 
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C.  Application  of  ARMAX  Method  to  Actuator  Number  Three 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX 
parameter  estimation  method  in  modeling  actuator  number  three.  The  ARX  based 
structure  [na=7  nb=6  nc=2  nk=l]  was  used 
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Zero  Pole  Plot  ARMAX  [7  6  2  1]  Model  *Data  Prefiltered* 
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Figure  C.6 
Zero-Pole  Plot 
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D.  Application  of  ARMAX  Method  to  Actuator  Number  Four. 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX 
parameter  estimation  method  in  modeling  actuator  number  four.  The  structure  [na=7  nb=6 
nc=2  nk=l]  was  used 
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Frequency  Response  Comparison 
Plot  of  ARMAX  MODEL[7  6  2  1]  Output  vs.  Actual  output 
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Residual  of  ARMAX  MODEL[7  6  2  1] 
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Zero  Pole  Plot  ARMAX  MODEL[7  6  2  1] 
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Zero-Pole  Plot 
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E.  Application  of  ARMAX  Method  to  Actuator  Number  Five 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX 
parameter  estimation  method  in  modeling  actuator  number  five.  The  structure  [na=7  nb=6 
nc=2  nk=l]  was  used. 
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Plot  of  ARMAX  MODEL[7  6  2  1]  Output  vs.  Actual  output  *Data  Prefiltered' 
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Zero  Pole  Plot  ARMAX  [7  6  2  1]  Model  *Data  PrefilterecT 
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F.  Application  of  ARMAX  Method  to  Actuator  Number  Six. 

The  following  plots  were  used  in  evaluating  the  application  of  the  ARMAX 
parameter  estimation  method  in  modeling  actuator  number  six.  The  structure[na=7  nb=6 
nc=2  nk=l]  was  used. 
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Zero  Pole  Plot  ARMAX  [7  6  2  1]  Model  *Data  PrefilterecT 
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